Monte Carlo Simulation: A Tool to Analyse
Magnetic Properties

Joan
Cano
and
Yves
Journaux

6.1 Introduction
The quest for molecular based magnets [1–3] and high-spin molecules [4–6], in
the wider context of molecular crystal engineering, has led to the synthesis of aesthetic extended networks [7] and high-nuclearity metal complexes [8, 9]. These
compounds give rise to interesting magnetic properties, such as spontaneous magnetization, or slow magnetic relaxation times and quantum tunnelling phenomena
[10–12]. Furthermore, the large majority of compounds belonging to these families
of materials often crystallise in novel topologies. In order to establish a correlation
between the structure and magnetic behavior of the compounds, it is essential to
develop suitable models for the description of the low lying and excited spin energy
levels. Unfortunately, the huge (or infinite) number of possible configurations in
these systems precludes the calculation of the exact partition function. As a consequence,thederivationofimportantthermodynamicpropertiessuchasthemagnetic
susceptibility and specific heat capacity cannot be done. This situation is typical
of systems studied by statistical physics which deals with systems with many degrees of freedom. Exact analytical theories are available in rare cases and in order
to tackle the calculation of thermodynamic properties, physicists have developed
approximate methods such as high temperature expansion of the partition function
[13], closed chain computational procedure [14, 15] or density matrix renormalizationgroupapproach(DMRG)[16,17].However,alltheseapproachesareoflimited
application or lead to uncontrolled errors which make improvement of the accuracy
of the results difficult. Monte Carlo simulation is the obvious choice to overcome
these problems [18]. The sources of errors are well known and the accuracy of the
calculation can be increased, in principle, by using more sample configurations and
by expanding the size of the simulated systems [18]. Furthermore, this approach
can be used for systems where analytic methods do not work. However, although
the Monte Carlo approach can be applied to many magnetic systems with different
types of interactions between the magnetic centers, this method remains simple to
program and affordable in term of computer power only in the case of the Ising
model [19] and the classical spin approximation (
=∞) [20]. Recent examples

Magnetism:
Molecules
to
Materials
V. Edited by J.S. Miller and M. Drillon
Copyright
c2005 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim
ISBN: 3-527-30665-X


190 6 Monte Carlo Simulation: A Tool to Analyse Magnetic Properties

in the literature provide interesting systems where these models can be applied.
In order to illustrate the power and the efficiency of Monte Carlo simulation in
molecular magnetism we will reduce the scope of this chapter to isotropic systems
with spin
=5/2 for which the classical spin approximation is satisfactory.

6.2 Monte Carlo Method
6.2.1 Generalities
Aclassicalprobleminstatisticalphysicsisthecomputationofaveragemacroscopic
observables such as magnetization
for a magnetic system. In the canonical
ensemble, the average magnetization
is defined as:

Mie.Ei/kT

M
i=1
(6.1)

.Ei/kT

e

i=1

It is generally not possible to compute exactly this quantity in Eq. (6.1) due to the
mathematical difficulties and the infinite or huge number of configurations. The
basic idea of Monte Carlo simulation is to get an approximation of Eq. (6.1) by
replacing the sum over all states with a partial sum on a subset of characteristic
states:

Mie.Ei/kT

M
i=1 (6.2)

.Ei/kT

e

i=1

In the limit, as
→∞, the sum formula of Eq. (6.2) equates to Eq. (6.1).
The first possible approach involved the random selection of the states for the
subset,i.e.adoptionofthesimplesamplingvariantofMonteCarlosimulation.This
approachhowever,hasmajordrawbacksastherapidlyvaryingexponentialfunction
in the Boltzmann distribution causes most of the chosen states to bring a negligible
contribution to Eq. (6.2). In order to get sensible results, the ideal situation would
be to sample the states with a probability given by their Boltzmann weight. As
will be shown below, this can be done by using the importance sampling approach.
Comparison between simple and importance sampling can be illustrated by the
fictitious system of 40 independent particles allowed to occupy 100 levels equally


6.2 Monte Carlo Method 191
100
80
60
40Energy20
Fig. 6.1. Occupation at 10 K of 100 energy levels equally spaced by 1 K using
0 a random selection (horizontal bars),
an importance sampling (dots) and a
population Boltzmann’s distribution (line).

spaced by 1 K. In Figure 6.1 is depicted the repartition of the independent particles
among the 100 levels at 10 K using a random selection and importance sampling
approach (Metropolis algorithm [21]). These two repartitions are compared to the
Boltzmann distribution.

This plot clearly shows that the high energy particles are too numerous in the
random sample when compared to the ideal Boltzmann distribution and will bias
the calculation of the average quantities. This is not the case for the sample obtained with the Metropolis algorithm. Even with a small number of configurations

(900) the repartition in the average sample is very close to the ideal Boltzmann
repartition. For a large number of configurations (900,000) the repartition of the
40 particles obtained by the Metropolis approach is indistinguishable from the
Boltzmann repartition.
The calculated average energies are 10.15 and 10.58 K for the random and the
Metropolis samples respectively (E=10.50 K for the 900,000 configurations
sample). The average energy calculated with the simple sampling is a poor approximation to the real average energy E=10.50 K. On the other hand, the
importance sampling approach gives sensible results, therefore it seems essential
to use this sampling method in Monte Carlo simulation [18]. In this approach,
the calculation of the average physical quantities is done by a simple arithmetic
average (Eq. (6.3))

0 1 2 3 4
T=10 K
1

M
Mi
(6.3)

i=1

But the configurations used in the arithmetic average are chosen according to
their Boltzmann weights. That is, for low temperature there are more low energy
configurations than high energy ones. Although the method looks reasonable, it
seems difficult to calculate the sampling probability p(Ci)
of a configuration Ci
which depends on the partition function (ZN
, Eq. (6.5)), that we are unable to
calculate

.Ei/kT

e

p(Ci) (6.4)

ZN


192 6 Monte Carlo Simulation: A Tool to Analyse Magnetic Properties

ZN e.Ei/kT
(6.5)

i=1

This tour de force can be accomplished by applying the Metropolis algorithm [21].

6.2.2 Metropolis Algorithm
The idea advanced by Metropolis et al. [21] is to generate each new configuration
Cj
from the previous one Ci
and to construct a so-called Markov chain [18]. The
probability of getting Cj
from Ci
is given by a transition probability W(Ci Cj).
It is possible to relate this transition probability W(Ci Cj)
to the probability of
a configuration p(Ci)
byconsideringthedynamicsoftheprocess.Atthebeginning
of the process, the probability of a configuration depends on the computer time
(number of iteration). Therefore, it is possible to calculate the probability of a
configuration Ci
at the time
1 through the following relation

p(Ci,
1) p(Ci,t)

W(Cj Ci)p(Cj,t) W(Ci Cj)p(Ci,t)
(6.6)
Cj
=Cj

After several iterations (thermalization process) the probability of a configuration
p(Ci,t)
must be independent of the computer time, that is

p(Ci,
1) p(Ci,t)
(6.7)

One possibility to cancel the second term of Eq. (6.6) is the so-called detailed
balance condition

W(Cj Ci)p(Cj,t) W(Ci Cj)p(Ci,t)
(6.8)

which can be rewritten in the case of the Boltzmann distribution for p(Ci)
as

W(Ci Cj)
p(Cj)
e.E(Cj)/kT
/ZN
e.E(Cj)/kT
W(Cj Ci) p(Ci) e.E(Ci)/kT
/ZN e.E(Ci)/kT
(6.9)

This condition gives a relationship between the ratio of the transition probabilities
and the ratio of the configuration probabilities. It is worth noting that Eq. (6.9) is
independent of the partition function ZN
and that all the quantities in the last ratio
of Eq. (6.9) are known or can be calculated. The next step is to give arbitrary values
to W(Ci Cj)
and W(Cj Ci)
respecting the detailed balance condition. In
1953, Metropolis, Teller and Rosenbluth proposed the simple following choice for
[21]

W(Ci Cj) e.E/kT
if 
0
(6.10)
1.0 if 
0

with 
E(Cj) E(Ci)


6.2 Monte Carlo Method 193
This choice satisfies the detailed balance condition and, more importantly, it can
be shown by simple arguments that a sequence of configurations generated by this
procedure represents a configuration sample according to the Boltzmann distribution [18]. Finally, the last step in a Monte Carlo simulation is to define whether the
new configuration is accepted to calculate average quantities from Eq. (6.3). According to the Metropolis algorithm, only the probability of the transition to a new
configurationisgiven,butnomoredirectinformationonthisconditionisprovided.
So, the success of a transition is ruled by a comparison of its probability with a real
random number
uniformly distributed between zero and unity (
∈[0,1]). Thus,
only when W(Ci
→Cj) is the new configuration accepted. This option is
sensible since most of the high-energy configurations will be rejected, especially at
low temperatures, where the transition probability W(Ci
→Cj)
reaches smaller
values than most random numbers r.

Although Monte Carlo simulation using the Metropolis algorithm appears to
be a simple alternative for the calculation of average quantities, some points are
delicate and can lead to unreliable results. The main points to be checked in order
to obtain a robust simulation are: the thermalization process, the size of the model,
the number of MC iterations and the random number generators.

6.2.3 Thermalization Process
Before calculating a physical observable, it is necessary to check that the memory
of the initial state is lost and the equilibrium distribution is reached, that is, the
probability of a configuration must be independent of the “computer time” (the
number of Monte Carlo steps, MCS) and should only depend on its energy. The
necessary time to get closer to the equilibrium can be very large at temperatures
lower than that of the magnetic ordering temperature. Sometimes, 3 ×104 MCS
site.1 are not enough to reach equilibrium. When a sample is in equilibrium at
highertemperatures(300K)andsuddenlycooled,theinitialconfigurationisfrozen
and a very large time is required to reach equilibrium in the new conditions. Shorter
times are required when a gradual decrease in temperature occurs. For example, in
a 3D cubic lattice the equilibrium is not completely reached after 105 MCS site.1
at 0.1 K, see Figure 6.2.

Toavoidthisproblemofslowrelaxationtowardequilibrium,twokeypointsmust
be considered. First, at each temperature, the configurations found at the beginning
of the simulation (first Monte Carlo loops) must be excluded in the calculation of
the physical observable. Generally, we discard the first 10% of configurations
generated by the MC algorithm, where equilibrium has not been reached. Second,
starting from a high temperature, a low cooling rate must be chosen according to
the following equation:

Ti+1 =kTi,
with 0.9 ≤k<
1.0 (6.11)


194 6 Monte Carlo Simulation: A Tool to Analyse Magnetic Properties
0
2
4
6
0 2 104 4 104 6 104 8 104 1 105
2 K
1 K
0.25 K
0.1 K
( T)
cal
/ ( T)
eq
MCS / site
Fig. 6.2. Limit number of MCS to reach
equilibrium as a function of the temperature
in a 3D cubic lattice when it is suddenly
cooled.
From this equation, the points at low temperatures, where the relaxation time is
large, get closer. Therefore, in the last Monte Carlo steps, when the sample has
reached equilibrium, a configuration is chosen as the initial configuration for the
next temperature. So, this configuration is placed close to the equilibrium condition.
6.2.4 Size of Model and Periodic Boundary Conditions
Except in the case of high nuclearity complexes, which have a finite size, it is
not possible to simulate real networks. Since the time of calculation is infinite in
the last cases, finite models must be considered in order to study the extended
network. On the other hand, it is necessary to use systems large enough to avoid
finite size or border effects [18]. To illustrate this point, let us take the example of
the antiferromagnetic S = 5/2 regular chain, where the exact law for a classical
spin approach is known [20]. The results of the MC simulation of χ|J | (magnetic
susceptibility) as a function of T/|J | for an increasing number of spins in the chain
showthat below100 spins the simulations are not accurate enough, so it is necessary
to reach 200 spins to avoid boundary effects at low temperatures (Figure 6.3).
0 20 40 60 80 100 120
4
12
20
40
100
400
Fisher
0.025
0.030
0.035
0.040
M
/ cm3mol-1
T / K
Fig. 6.3. χ|J | versus T/|J | plot for a series of linear systems with an increasing number of
sites. The results are compared with Fisher’s law for a one-dimensional system of classical
spin moments [20].
χ χ χ

6.2 Monte Carlo Method 195
Fig. 6.4. Illustration of how periodic
boundary conditions are used to diminish
the model size without introducing border
effects.
Although it is possible to simulate a chain of 200 spin moments within a reasonable
time, itwould take too much time to simulate a 3Dnetwork on a 200×200×200
model. In fact, the threshold size to avoid finite size effects over a wide range for a
3D system is smaller, but the required size is still too large to allow an affordable
calculation time. The periodic boundary conditions (PBC) are used to diminish the
model size without introducing border effects [18]. Thus, for instance, the first and
last spin moments of a chain are considered as nearest-neighbours and, in consequence,
all spin moments become equivalents, see Figure 6.4. These PBC can be
extended to 2d and 3d networks as shown below.
In a one-dimensional system, the PBC conditions reduce considerably the finite
size effects, so 20 spin moments are enough to obtain a nearly perfect simulation
of the magnetic behavior (Figure 6.5).
0.020
0.025
0.030
0.035
0 20 40 60 80 100 120
4 cyclic
8 cyclic
20 cyclic
40 cyclic
Fisher
..M / cm3mol-1
T/K
Fig. 6.5. χ|J | versus T/|J | plot for a series of cycles
with an increasing number
of sites. The results are compared
with Fisher’s law for
a one-dimensional system of
classical spin moments [20].

196 6 Monte Carlo Simulation: A Tool to Analyse Magnetic Properties

In practice, the model size considered in the Monte Carlo simulations is double
the minimum size for which border effects are absent.

6.2.5 Random Number Generators
The use of random numbers is the core of Monte Carlo simulations. Thus, finding
a good random number generator (RNG) is a major problem. In an ideal situation the random numbers would be generated by a random physical process. In
practice, computers are used to carry out this function using mathematical subroutines. Actually, these generated numbers are not random and are referred to
as pseudorandom numbers. However, this difference is not very meaningful if a
RNG satisfies some important criteria. In general, RNGs supplied with compiler
packages are dirty generators, so it is necessary to find a proper subroutine. A good
introduction to RNGs is given in Ref. [22]. All mathematical RNGs supply a finite
number sequence, which must be reproducible in any computer. The period of this
sequence must be long and, at least, very much larger than the required random
numbers sequence to simulate a physical property at a certain temperature. On the
other hand, the number produced by a RNG must be apparently random, in other
words,theymustpresentahomogeneousdistribution,avoidingnumbersturningup
as a series of sequences involving numbers of a similar magnitude. To summarize,
these RNGs must satisfy the statistical tests for randomness. Unfortunately, some
RNGs fulfil statistical tests but fail on real problems. Thus, it is necessary to check
the RNGs on real problems that have already been solved. Some authors suggest
testing the RNGs using the MC methods in the calculation of energy for a 2D Ising
network [23].

6.2.6 Magnetic Models
Althoughthenatureoftheinteractionbetweenthemagneticionsiselectrostatic,the
magnetic data can be well described using effective spin hamiltonians reminding
of a magnetic interaction. Theoreticians have justified the use of such hamiltonians
for magnetic systems. Most of the studies have been based on the Heisenberg and
Ising Hamiltonians, which can be written in general as:

=
Jij
Siz.Sjz
+Jij
⊥(Six.Sjx
+Siy.Sjy)
(6.12)

j>

where Sik
are components of the spin vectors Si
, and Jij and Jij are the exchange
coupling constants. The Ising and Heisenberg models correspond to cases where
Jij
⊥
0 and Jij

Jij
⊥, respectively. The Ising model is adapted to strongly
anisotropic ions, but in spite of its mathematical simplicity nobody has been able


6.2 Monte Carlo Method 197
to solve it exactly beyond the 2D square lattice [19]. On the other hand, the Heisenberg model is adapted to isotropic systems, but it is not possible to solve it except
for some finite systems. However, Monte Carlo simulation is a useful tool to describe the magnetic properties of systems where exact solutions are not known. We
have shown that the Metropolis algorithm allows one to sample the configurations
according to the Boltzmann distribution. The core of the algorithm is the comparison of a random number with the quantity e.E/kT
. In theory, it is necessary to
diagonalize the full energy matrix built from the Heisenberg hamiltonian to know
the energies of the configurations that allow the calculation of e.E/kT
. Thus, apparently, it seems that we have gone back to the starting point. It is possible to
overcome this problem by using a Quantum Monte Carlo approach but this is beyond the scope of this chapter [24–27]. There are many interesting compounds that
containionswithspins
≥2 (Mn(II) or Fe(III)), where there is another possibility.
It has been shown that these spins can be considered as classical vectors. However, in order to compare the calculated values with experimental observations, the
classical spin vectors are scaled according to the following factor:

Si

Si(Si
+1)
(6.13)

With this approximation, the Heisenberg hamiltonian is reduced to

=
Jij Si(Si
+1) Sj(Sj
+1)
·cos θij
(6.14)
j>

which allows one to easily calculate the configuration energies requested for the
Metropolis algorithm (CSMC method). Thus, this chapter focuses on the
=5/2
systems, where the classical spin approach can be used.

6.2.7 Structure of a Monte Carlo Program
All the ingredients to write a Monte Carlo program are available. An abstract of
this program is shown in Figure 6.6.

Two remarks must be made: (a) the initial spin configuration of the network or
cluster is chosen randomly, but other choices are possible; and (b) the sites of the
network are not explored randomly for the spin orientation update but systematically through a loop. It has been shown that this approach gives good results for
equilibrium configurations.

After the generation of the sample using the Metropolis algorithm all the thermodynamic quantities can be calculated. It has been shown that the magnetization
is calculated as the simple arithmetic average

1

M
Mi
(6.15)

i=1

And it is also possible to calculate the average energy:


198 6 Monte Carlo Simulation: A Tool to Analyse Magnetic Properties

Initialisation of the network
Random orientation for the spins

T=T initial

..E ≤0
Tnew = T xK
Calculate ..E
..E > 0
random number R
R ≤ e-..E/kT
keep the old spin orientationtake the new spin orientation
i = maxsite
istep >
thermalisation steps
memorize configuration
istep = maxstep
calculate average
quantities .... M
i = i +1
istep = 0
site i= 0
Yes No
istep =istep+1
The new spin components
Sx(i) , Sy(i), Sz(i) are choosed
randomly
Yes No
No Yes
Yes No
T =Tend
Yes
No
STOP
Fig. 6.6. Flow chart of a Monte carlo program using the classical spin approach.

1

E
Ei
(6.16)

i=1

as well as the magnetic susceptibility and specific heat, which are calculated as the
fluctuation of magnetization and energy, respectively.

χ=M2.M2;Cp=E2.E2 (6.17)


6.3 Regular Infinite Networks 199
6.3 Regular Infinite Networks
In order to fully understand and fine-tune the physical properties of magnetic
materials, it is necessary to gain as much information as possible, such as the g-
factors or the interaction parameter
between the magnetic ions. For a simple
system, it is possible to get the values of these parameters by fitting a theoretical
model to the experimental data. So, the calculation of the magnetic susceptibility
is often combined with a least-square routine allowing the determination of the
best parameters. In practice, a least-square fit by Monte Carlo simulation takes a
lot of computer time. Nevertheless, for networks with only one or two interaction
parameters, empirical laws using reduced variables can be established from Monte
Carlo simulations [28]. The magnetic susceptibility can be given by an expansion
function

aJ

fb
+ε(H)
(6.18)

TT

When the magnetic field is close to zero the ε(H)
term is negligible and the magnetic susceptibility becomes field independent. In this case, there is only one χ|J|versus T/|J|curve for all
values. So, it is possible to obtain empirical laws from
the Monte Carlo simulations which depend on the reduced temperature
=T/|J|.
These empirical laws, which have been derived for several regular networks (1D,
square and honeycomb 2D and cubic 3D), take the form:

a0
aiβi

2

χ|J|=
gi=1 (6.19)

k+1

4
1
bjβj

j=1

The coefficients associated with the highest degree of the polynomials for both
the denominator and the numerator are set so that they converge to the Curie law
at high temperatures (χT 4.375 cm3 K mol.1, for
2). Furthermore, the
zero-grade terms in the numerator are fixed so that they converge to the finite
χ|J|values obtained by the simulations at low temperature. The exact numerical
coefficients associated with the empirical laws derived for cubic, diamond and 3connected 10-gon (10, 3) 3D networks are given in Table 6.1, and those for square
and honeycomb 2D networks are shown in Table 6.2. Empirical laws can also be
found for alternating systems with different magnetic interactions, but they present
a more complicated form. Equations for other systems that are not shown in the
present manuscript are available from the authors.

Acomparisonismadewiththehightemperatureseriesexpansionofthepartition
function (HTE) for 2D honeycomb and square and 3D cubic networks by Stanley


200 6 Monte Carlo Simulation: A Tool to Analyse Magnetic Properties

Table 6.1. Coefficients of the rational functions providing the thermal variation of the reduced
magneticsusceptibility χ|J|
as a function of
T/|J|
for simple cubic, diamond and (10, 3)
cubic networks (
in K) (Eq. (6.19)) [36].

Coefficient Cubic Diamond (10,3) Cubic

a0 0.0815865 0.116 0.156

b0111
a1000
b1000

a21.22599
10.5 1.85958 1.20777
10.4

b2 .2.78782
10.3 11.7921 5.11803
10.5
a3000
b3000

a4 .5.34657
10.7 324.872 6.12417
10.4

b49.71169
10.6 2795.33 .2.47313
10.4
a5000
b5000

a63.56382
10.8 2.5012 .1.0435
10.6

b61.37954
10.9 18.033 .2.32108
10.4
a7000
b78.14585
10.90 0

a8 0 0.264794 3.99986
106

b8 0 0.816704 8.05117
10.6
a9000
b9 0 0.0605244 9.14253
10.7

et al., and Lines et al. and with Fisher’s law for a chain (Figure 6.7) [13, 20, 29–
33]. Agreement between MC simulations and the other approaches is excellent at
temperatures higher than that of the maximum value of χ|J|. Nevertheless, below
this temperature there is a noticeable discrepancy since the HTE method is not
applicable in this region, whereas there is a perfect agreement between the MC
simulation and Fisher’s law for a regular chain over the whole temperature range.
On the other hand, as expected, the maximum value of χ|J|
increases and its
position is displaced towards lower temperatures when the dimensionality of the
networkandtheconnectivitybetweenthemagneticionsdecrease,sincethenumber
of spin correlation paths also decrease [28, 34]. It must be noticed that for the 3D
system the maximum corresponds to the antiferromagnetic ordering temperature,
whereas for the 1D and 2D networks it is well established that there is no magnetic
ordering for the Heisenberg model.

We have tested the CSMC approach to fit the magnetic data for
[N(CH3)4][Mn(N3)3] [35, 36], which crystallizes in a regular cubic network


6.3 Regular Infinite Networks 201
Table 6.2. Coefficients of the rational functions providing the thermal variation of the reduced
magnetic susceptibility χ|J|as a function of
T/|J|for square and honeycomb 2D
networks (
in K) (Eq. (6.19)).

Coefficient Square Honeycomb
a0 .121201.0 2.82178
b0 .1.05473 ×106 .582.803
a1 311085.0 82.6317
b1 2.72275 ×106 3830.22
a2 .289512.0 .110.786
b2 .2.56424 ×106 .8491.63
a3 .117474.0 .248.245
b3 1.07481 ×106 8711.63
a4 .19202.0 626.252
b4 .201403.0 .4401.53
a5 358.413 .530.795
b5 17648.1 982.435
a6 428.864 220.820
b6 .275.529 .34.5723
a7 .73.9798 .47.0672
b7 .39.9761 .12.2174
a8 4.375 4.375
b8 .5.4875 .2.11426
a9 0.000 0.000
b9 1.000 1.000

(Figure 6.8). Its magnetic behavior can be reproduced using
=.5.2cm.1 and
2.025 [36]. These values are close to those found with the HTE method
(
=.5cm.1) [35]. It is worth noting that the agreement between the MC simulation and the experimental data is very good, even at low temperatures, confirming
the classical behavior for this 3D network.

An interesting comparison between the magnetic behavior of three different
antiferromagnetic regular 3d networks is shown in Figure 6.9. These 3D systems
correspond to primitive cubic, diamond and 3-connected 10-gon (10, 3) cubic
networks. As expected, the antiferromagnetic ordering temperature TN/|J|is displaced toward a lower temperature as the connectivity between the magnetic sites
decreases [36]. Below TN/|J|, in an ordered phase, our results could be compared
to the less accurate mean field approximation. In this approach, as in our MC simulations, the expected limit of the χ|J|value at T/|J|=0 is equal to 2/3 of its
maximum value. So the CSMC method is able to reproduce the physical behavior
in the paramagnetic and in the ordered phases, while the mean field approximation


202 6 Monte Carlo Simulation: A Tool to Analyse Magnetic Properties
(a) (b)
(c) (d)
0.050
0.100
0.150
0.200
0.250
0.300
0 10 20 30 40 50
..|J| / cm3mol-1
T / |J|
(a)
(c)
(b)
(d)
Fig. 6.7. χ|J | versus T/|J | plots obtained by MC simulation for 1D, 2D honeycomb, 2D
square and 3D cubic networks in a Heisenberg model. These plots are compared with those
obtained by the high temperature expansion method and Fisher’s law [13, 20, 29–33].
0.010
0.015
0.020
0.025
0.030
0.035
0 50 100 150 200 250 300
.. / cm3mol-1
T / K
Fig. 6.8. Crystal structure and magnetic properties of [N(CH3)4][Mn(N3)] [35, 36]. The experimental
data (circles) and the simulations by Monte Carlo methods (solid line) are shown.
leads to a large overestimation of the ordering temperature and the HTE method is
limited to the paramagnetic region [34].
In [FeII(bipy)3][MnII
2 (ox)3], a compound previously described by Decurtins et al.
[37], where bipy=2, 2-bipyridine and ox=oxalate, the Mn(II) ions are connected
via oxalate bridging ligands to build up a three-dimensional 3-connected 10-gon
network (Figure 6.10). From CSMC simulations, an antiferromagnetic interaction
is found for this compound with J = .2.01 cm.1 [36]. This value is in agreement
with those found in the literature for other dinuclear complexes and regular chains
incorporating oxalate groups as bridging ligands.

6.4 Alternating Chains 203
0.050
0.10
0.15
0.20
0.25
0 10 20 30 40 50
..
M
|J| /cm3 mol -1 K
T / |J|
(c)
(b)
(a)
(a) (b)
(c)
Fig. 6.9. MC simulations of χ|J|
versus T/|J|
plots for three different antiferromagnetic
regular 3D networks: (a) primitive cubic, (b) diamond and (c) 3-connected 10-gon (10, 3)
cubic networks [36]. The line without symbols represents the theoretical curve found by the
HTE method.

0.022

/ cm3mol-1

M

0.018

0.014

0.010
100 200 300
T / K

Fig. 6.10. Crystal structure and magnetic properties of [FeII(bipy)3][MnII2 (ox)3] [37]. The
experimental data are shown as circles and the Monte Carlo and HTE simulation as bold and
normal lines, respectively [36].

6.4 Alternating Chains
Alternating
5/2chainswithtwoormoredifferentexchangecouplingconstants
have also been investigated. An interesting example is that of a chain presenting
an interaction topology J1J2, that is, two different consecutive interactions (J1 and
J2) that repeat along the chain (...J1J2J1J2J1J2...). Drillon et al. have derived
an exact analytical law in the frame of the classical spin approach to analyse the
magnetic behavior of these systems [38]. On the other hand, the versatility in the
coordinationoftheazidoligandledtoseveralinteractiontopologies.Inthe[MnII(2pyOH)2(N3)2]
compound (2-pyOH
2-hydroxypyridine) the manganese(II) ions
are connected by μ-1,3-azido bridging ligands (Figure 6.11) [39]. It is well known


204 6 Monte Carlo Simulation: A Tool to Analyse Magnetic Properties

Fig. 6.11. Crystal structure and magnetic properties of [MnII(2-pyOH)2(N3)2]
[39]. The experimental data are shown as circles and the simulations by Monte Carlo method and Drillon’s
law as solid and dashed lines, respectively [38].

Fig. 6.12. Crystal structure and magnetic properties of [MnII(bipy)(N3)2]
[40, 41]. The experimental data are shown as circles and the simulations by Monte Carlo method and Drillon’s
law as solid and dashed lines, respectively [38].

that this kind of bridge leads to antiferromagnetic interactions. However, a μ-1,1azido bridging ligand is also present in the [MnII(bipy)(N3)2]
compound; so, a
ferromagnetic interaction is expected in this case (Figure 6.12) [40, 41].

The experimental data were simulated by the CSMC method and Drillon’s law,
and a good agreement is found between both methods and the experimental data
(Table 6.3).

The compound [Mn(Menic)(N3)2]
(Menic
methylisonicotinate) represents a
more complex alternating 1D system. In this chain, a 1:4 ratio for the μ-1,3-and μ1,1-azidobridgingligandsconnectingthemanganese(II)ionsisfound(Figure6.13)
[42]. Nevertheless, a more complicated interaction topology than J2J2J2J2J1 is
observed for this compound. In this way, as there are two different MnNazidoMn
bond angle (α)
values for the [Mn2(μ-1,1-azido)2] entities (101.1o
and 100.6o
,
Figure 6.13), a J2J3J3J2J1 interaction topology must be considered.

In a recent paper, Drillon et al. conclude that several alternating ferroantiferromagnetic homometallic one-dimensional systems present similar magnetic behavior to that of the ferrimagnetic chains [43], which are described by


6.4 Alternating Chains 205
Table 6.3. Best parameters obtained by fitting a theoretical model to the experimental data
for [MnII(2-pyOH)2(N3)2]
and [MnII(bipy)(N3)2]n. The fits have been performed using the
CSMC method and Drillon’s law [38–41].

Compound Method
J1/cm.1 J2/cm.1
[MnII(2-pyOH)2(N3)2]
MnII(bipy)(N3)2]
MC
Drillon
MC
Drillon
2.04
2.03
1.98
1.99
.13.2
.13.8
.12.9
.12.9
.12.3
.11.7
+4.9
+5.0

Kahn as systems that contain two different near-neighbour spin moments antiferromagnetically coupled [15].

[Mn(Menic)(N3)2]
constitutes a beautiful example of these systems. Thus, its
χT
versus
experimentalcurvepresentsaminimum.Moreover,amaximumisalso
observed at lower temperature, which is characteristic of this particular interaction
topology.Onlythe J2J3J3J2J1 model provides a correct description of the magnetic
behavior, even at low temperatures (Figure 6.13). A good fit is obtained with the
set of parameters J1 =.15.6cm.1, J2 =1.06cm.1 and J3 =1.56cm.1. These
results agree perfectly with those obtained from a proposed exact analytical law
[42]. As is known, the
parameter and the
angle are related. Thus, the J2 and J3
valuesareinagreementwiththetheoreticalmagneto-structuralcorrelationfoundby
Ruiz et al. [44], supporting the consideration that two different exchange coupling
constants for the [Mn2(μ-1,1-azido)2] entities has a physical meaning and is not
the result of a mathematical artifact.

100.6o 101.1o
J2 J2J3 J3 J1
Fig. 6.13. Crystal structure and magnetic properties of [Mn(Menic)(N3)2]n. The experimental
data are shown as circles and the simulations by Monte Carlo method and the exact analytical
law proposed by us as solid and dashed lines, respectively. The interaction topology is shown
in the picture of the crystal structure [42].

The presence of a minimum and a maximum in the χMT
versus
curve, can be
explained by considering instant spin configurations at several temperatures provided by the MC simulation process (Figure 6.14) [42]. In this way, the stronger
antiferromagnetic coupling promotes an antiparallel spin configuration and the


206 6 Monte Carlo Simulation: A Tool to Analyse Magnetic Properties

Fig.6.14. Thermalvariationofthespinconfigurationforanalternating J2J3J3J2J1 interaction
topology. The antiferro-and ferromagnetic interactions are represented by bold and dashed
lines, respectively.

χMT
product decreases on cooling. At lower temperatures, in spite of the weak
character of the ferromagnetic interactions, as they are present in a greater proportion (4:1), the ferromagnetic alignment of the resulting spins becomes efficient.
Finally, the strongest antiferromagnetic interaction dominates the magnetic behavior, and χMT
increases to reach a maximum then further decreases attaining a zero
value at 0 K.

6.5 Finite Systems
In a system where the local spin moments have an infinite value (Si
=∞) the
number of microstates (Sz
values) is infinite, as in a classical spin approach where
the spin vector can be placed along infinite directions. Thus, the classical spin
approach is less correct when the value of the local spin moment decreases, since
the quantum effects become non-negligible. We have verified previously that the
classical spin approach can be used to analyse the magnetic behavior of periodic
systems. In a real non-periodic system the number of states is limited, especially
for small systems with a few paramagnetic centers, and the energy spectrum is
far from being a continuum, so the mentioned quantum effects could be more
important in these systems. The question is whether the magnetic behavior of these
non-periodic systems can be reproduced by MC simulation within the framework
of the classical spin approach. In other words, does the applicability of the classical
spin approach depend only on the values of the local spin moments or does the size
of the network have some influence? Moreover, is it possible to simulate a discrete


6.5 Finite Systems 207
(a) (b)
J 0.500

1

0.400

J

2

0.300

J

0.200

3

..|J| / cm3mol-1K

0 3 6 9 12 15
..
..
..
3
2
1
T / J

Fig. 6.15. A comparison between the theoretical χ|J|
versus T/|J|
plots simulated from the
exact quantum solution (symbols) and Monte Carlo methods (lines) for a series of small linear
models.

spectrum from a continuous energy spectrum? For any system at low temperatures
or for very small systems this task can be especially difficult since there are few
populated states. In this section, the limits where the classical spin approach can be
applied will be established. Thus, from the study of some systems where an exact
quantum solution is available (Figure 6.15), it is possible to check the validity of
the classical spin approach. In Figure 6.15, the curves for several linear systems
obtained by CSMC simulation are compared to those calculated from an exact
quantum method. The classical spin approach is not valid at low T/|J|
values, due
to the small number of populated states, which can be considered as a quantum
effect. Also, from Figure 6.15 it can be concluded that the higher the number
of paramagnetic centers, the lower the quantum effect. Thus, for any system at
T/|J| 4, the classical approach can be applied, and it is valid for a wider
range of T/|J|
as the number of paramagnetic centers increases. So, for more
extended systems where the exact numerical solutions cannot be calculated, the
MC simulation in a classical spin approach will be a powerful tool to study their
magnetic behavior.

The same comparison between MC simulations and exact quantum numerical
solutions has been made for spin topologies presenting more than one coupling
constant. Two examples are shown in Figures 6.16 and 6.17, where only interactions of antiferromagnetic nature are present. χT
versus T/|J|
curves have been
simulated for different J/
values. MC simulations are valid for T/|J|
values
higher than 1.5, for J>J, and similar conclusions as above are reached.

Moreover, several real complexes have been studied by a CSMC method. As
an example, in Figure 6.18, two interesting clusters containing ten and eighteen
iron(III) ions respectively, with a ring structure, named ferric wheels, are shown
[45, 46]. From the magnetic point of view, these clusters are beautiful examples
of systems that can be used as models for the interpretation of the magnetic prop


208 6 Monte Carlo Simulation: A Tool to Analyse Magnetic Properties
0.000
1.000
2.000
3.000
4.000
5.000
0 5 10 15 20
..T / cm3mol-1K
T / J
..
.. ..
..
.. ..
..
1.0
0.8
0.6
0.4
0.2
0.0
J
J'
Fig. 6.16. A comparison between the theoretical χT versus T/|J | plots as a function of the
J /J ratio (α) for the model shown in the picture. The plots have been simulated by exact
quantum solution (symbols) and Monte Carlo methods (lines).
0.000
2.000
4.000
6.000
8.000
0 5 10 15 20
..T / cm3mol-1K
T / J
..
.. ..
..
.. ..
..
1.0
0.8
0.6
0.4
0.2
0.0
J
J'
Fig. 6.17. A comparison between the theoretical χT versus T/|J | plots as a function of the
J /J ratio (α) for the model shown in the picture. The plots have been simulated by exact
quantum solution (symbols) and Monte Carlo methods (lines).
erties of linear chains. (versus T plots have been simulated by a CSMC method,
considering magnetic interactions only between nearest neighbours. The obtained
values of the coupling constants agree perfectly with those obtained by an exact
analytical classical spin law for 1D systems.
6.6 Exact Laws versus MC Simulations
In previous sections, it has been shown that the MC method is a high-performance
tool to simulate the magnetic behavior of many different systems. Notwithstanding,
in some cases exact classical spin (ECS) laws are also available, so experimental
data can more easily be processed. Thus, the question arises as to whether to use a

6.6 Exact Laws versus MC Simulations 209
g = 1.98
J = -9.8 cm-1
g = 1.985
J1 = -19.1 cm-1
J2 = -8.0 cm-1
…J1 J1 J2…
(a)
(b)
Fig. 6.18. Crystal structure, experimental (circles) and MC simulated (lines) magnetic properties of: (a) [Fe(OCH3)2(O2CCH2Cl)]10 and (b) [Fe(OH)(XDK)Fe2(OCH3)4(O2CCH3)2]6
(where XDK is the anion of m-xylylenediamine bis(Kemp’s triacid imide)) [45, 46].

MCmethodwhenananalyticallawcanbeapplied.Atthemoment,therearealready
ECS laws for several 2D networks. It is very important to understand how these
laws are elucidated and what are their applicability limits and this is the subject of
the present section. First, a method to obtain an ECS law for a 1D system, that is
the Fisher’s law [20], is described.

6.6.1
A Method to Obtain an ECS Law for a Regular 1D System:
Fisher’s Law
The evaluation of any physical property at a precise temperature requires the solution of two integrals (Eq. (6.1)): (a) the value of this property as the sum of the
contributions from each of the possible states, and (b) the normalization factor, that


210 6 Monte Carlo Simulation: A Tool to Analyse Magnetic Properties

θι
Fig. 6.19. Illustration of the angle between two coupled
i+1 vectors.

is given by the partition function (ZN
). Considering the spin moments as vectors,
the energy of any configuration is given by Eq. (6.14), where θij
is the angle between two coupled vectors. In the case of a regular chain, the partition function
can be described as

1
ZN sin θi
e(
cos θi)dθi
(6.20)
i=1 02

being equal to .JS(
+1)/
. It must be pointed out that there is a term, sin θ,
thattakesintoaccountthedifferentarrangementsthatcanbegeneratedataconstant
angle θ, that is, from a precession of the second vector referred to the direction of
the first one, see Figure 6.19.

The next equation is found from solving this integral


N

sinh(x)

ZN (6.21)

The magnetic susceptibility in zero field can be calculated from the total spin pair
correlation function, which can be defined as the sum of the individual spin pair
correlation functions:

2β2
NN

(
4kT sizsjz
(6.22)
i=0
=0

These functions, that provide the average arrangement of all spin moments
referred to one of them, can be evaluated in a similar way to the partition function.
Defining a pair correlation function by

3

si
·si+1
sizsiz
+1
(6.23)
ZN

In a 1D system the integrals may be factorised as before.

1

=si
·si+1
cos θi
sin θi
e(
cos θi)dθi
(6.24)

2
This factorisation involves an independent character for the spin pair correlation
function concerning only two near-neighbour centers. However, this is not the case

0


6.6 Exact Laws versus MC Simulations 211
for topologies other than a chain where this methodology cannot be so easily
applied. In this way, Langevin’s function is obtained, which describes how one
spin moment is placed with respect to its neighbours.

1

coth(x) (6.25)

Obviously, the spin correlation function of vector
with itself is unity. On the
contrary, the neighbouring vector
1 is correlated to vector
by Langevin’s
function (u). On the other hand, the neighbouring vector
2 is correlated to i
through vector i+1. Thus, the spin pair correlation function of vectors i+2and
is
2. From the summation in Eq. (6.22) and considering the obtained individual spin
pair correlations, the series shown in Eq. (6.26) is constructed. The factor 2 that
appears in some of the terms of the equation comes from the fact that an infinite
chain grows in the two directions of the chain axis. By expanding the summation
over integer
values, the wellknown Fisher’s law is obtained (Eq. (6.27).

n

χT χTfree–ion(1
2u 2u
2
2u
3
...
χTfree–ion 1
2u
(6.26)
n=1

1

χT χTfree–ion (6.27)

1

6.6.2 Small Molecules
The simplest case that can be studied is a system with only two paramagnetic
centers. Following the methodology detailed in the preceding section, the next
ECS law can be deduced

χT χTfree–ion(1
u)
(6.28)

The χ|J|
versus T/|J|
plots obtained from Eq. (6.28), from a CSMC simulation
and from the exact quantum solution are shown in Figure 6.20. A good agreement
between the three methods is found, and some discrepancies appear only at low
values of T/|J|. The MC method gives a better result than the ECS law for this
T/|J|
region, since an approach has been made in the calculation of the partition
function.

In this same way, several comparisons between the three methods have been
made for a series of similar models, and it can be concluded that there is a good
agreement amongst them. Nevertheless, there is no such agreement in systems
where the interaction topology presents closed cycles, as those shown in Figure 6.21. In some of these cases where the exact quantum solution is available,
it has been observed that ECS laws do not simulate the magnetic behavior of the
system properly. In the simplest case, that is, a triangle, vector 1 is correlated to


212 6 Monte Carlo Simulation: A Tool to Analyse Magnetic Properties
0.3
0.4
0.4
0.5
0.5
0.6
0 2 4 6 8 10
..|J| / cm3mol-1K
T / J
ECS law
CSMC
Quantum result
Fig. 6.20. Theoretical χ|J | versus T/|J | plots obtained by exact quantum solution, exact
classical law and CSMC simulation.
vector 3 either through the left or the right-hand ways, increasing the number of
correlation paths. Also, vector 1 can be correlated to itself through a correlation path
that involves a whole turn. This turn can be made clockwise or counter-clockwise,
but both paths are equivalent, so it must be considered only once in the ECS law,
as shown in the next equation.
χT = χTfree–ion(1 + 2u + 2u3 + u3) (6.29)
Notwithstanding, as vector 1 is correlated to itself through all other vectors, it is not
possible to split the spin correlation function in their individual spin pair correlation
functions. Thus, this methodology is not useful in these systems, and the derivation
of an ECS law is a hard and difficult task.
Series 1 Series 2
Fig. 6.21. Systems where the interaction topology presents closed cycles.

6.6 Exact Laws versus MC Simulations 213
a) b)

20.0

5.0

18.0

16.0

2 3 4 5 6 7 8
(T / J)
lim
(T / J)

lim

14.0

12.0

3.0 10.0
3 4 5 6
Cycle size
7 8
8.0
Number of cycles

Fig. 6.22. The limit value of T/
for perfect agreement between the exact classical and
quantum solutions increasing: (a) the cycle size and (b) the number of cycles (see Series 1 and
2 in Figure 6.21). The condition to control the quality of the agreement is stricter in case (b)
than in case (a) in order to facilitate the analysis of the results.

A study has been performed on the series of topologies shown in Figure 6.21.
From a comparison of the results of the ECS laws with those from MC simulations
it can be deduced that, by decreasing the number of triangular cycles (Series 2)
or increasing the size of the cycle (Series 1), the validity range of the ECS laws
increases and, consequently, the T/|J|
threshold decreases (Figure 6.22). In the
firstcase,thiseffectisduetoadecreaseinthenumberofcorrelationpathsinvolving
one or more closed cycles. In the second case, when the cycle size increases, the
value of the spin correlation function involving a closed cycle path becomes lower.
Therefore, this kind of correlation path is negligible in an infinite size ring, as a
result again obtaining Fisher’s law.

6.6.3 Extended Systems
Two analytical laws have been derived to date for a 2D network. These ECS laws
have been deduced by Cur′

ely et al. for alternating square and honeycomb networks
[47–50], where there is only one magnetic interaction along a chain and a different interchain interaction (Figure 6.23). From these equations, ECS laws can be
obtained for the corresponding regular networks. In this way, it could be expected
that more complex topologies or 3D networks could be solved, and that it would
not be necessary to use MC methods to simulate their magnetic behavior. On the
other hand, MC methods allow one to consider as many coupling constants and
g-factors as desired for any system, but this is not the only reason for continuing
to use MC methods.
From an analysis of Cur′

ely’s law for an alternating 2D square network it can
be observed that the spin correlation function is, surprisingly, the product of the
one-dimensional spin correlation functions along each spatial direction [50, 51]. In
this topology, when ferromagnetic and antiferromagnetic interactions are present it


214 6 Monte Carlo Simulation: A Tool to Analyse Magnetic Properties
J
J'
J
J'
2d square 2d honeycomb
Fig. 6.23. Alternating square and honeycomb networks.
is expected that χT will reach a zero value at 0K, whatever the magnitude of these
two interactions. Nevertheless, from Cur′ely’s law, in the case of a ferromagnetic
interaction stronger than the antiferromagnetic one, the χT product diverges on
cooling, and when both interactions are of the same magnitude, χT versus T
follows Curie’s law, that is, the system, surprisingly, behaves as if the spin moments
do not interact at all. Moreover, the coefficients of the high temperature expansion
for Cur′ely’s law do not agree with those obtained by Camp et al. or Lines [30,
33, 51]. These remarks, as some that will be made later, can also be extended
to Cur′ely’s law for 2D honeycomb networks. Thus, for instance, [Mn(ox)2(bpm)]n
presents an alternating honeycomb network, where the oxalate (ox) bridging ligand
acts as an exchange pathway along one of the directions and the bipyrimidine (bpm)
ligand connects the chains (Figure 6.24) [52]. Excellent fits of the model to the
experimental data for this compound have been obtained both fromMCsimulations
and from Cur′ely’s law [28, 49]. However, the J constant values obtained from MC
Fig. 6.24. Crystal structure and magnetic properties of [Mn(ox)2(bpm)]n [52]. The experimental
data (circles) and the simulations by Monte Carlo methods (solid line) are shown.

6.6 Exact Laws versus MC Simulations 215
simulations agree much better with those found in dinuclear and one-dimensional
systems with oxalate or bipyrimidine as bridging ligands [28, 52].
A more detailed analysis of the elucidation of an ECS law for a regular 2D
square network will allow one to find the limitations of this methodology. Thus,
as in this kind of network all paramagnetic centers are equivalent, only the total
spin correlation function referred to one spin moment must be evaluated. As has
been previously said, correlation of spin moment A to itself is 1, correlation of spin
moment B to A is given by the Langevin’s function (u), and correlation of spin
moment C to A is u3 (Figure 6.25). In this way, the summation in Eq. (6.30) is
generated, and its resolution leads to Cur′ely’s law.
χT = χTfree–ion..
1 + 2
∞
i=1
ui + 2
∞
j=1
uj + 4
∞
i=1
∞
j=1
ui+j..
= χTfree–ion 1 + u
1 . u2
(6.30)
Notwithstanding, three short correlation paths exist between the A and C spin
moments (1, 2 and 3), so the u3 term must be computed three times. An analytical
law can easily be deduced from the summation generated by computing all these
paths (even for the alternating case), but an infinite number of correlation paths
such as 4 and 5, and a finite number such as 6, have been omitted from this reasoning.
These longer correlation pathways are not so important for a wide range of
T/|J | values, but they are numerous and the contribution from all of them can be
significant and must not be disregarded. Therefore, the calculation of ECS laws for
more than 1D becomes impossible. In Table 6.4, the number of different spin cor-
A
B
C
1
2
A
C
C
5
3
4
A
C
C
6
Fig. 6.25. Illustration of correlation
paths to spin moments.

216 6 Monte Carlo Simulation: A Tool to Analyse Magnetic Properties
Table 6.4. Number of individual spin correlation paths (λ) as a function of the correlation path
length (n) in a 2D square network, obtained by SPPA method.
n λ n λ n λ
0 1 8 6674 16 22155058
1 4 9 18600 17 60555564
2 12 10 51480 28 165126324
3 36 11 142412 29 450294176
4 104 12 391956 20 1225587036
5 300 13 1078612
6 848 14 2956928
7 2392 15 8105796
0.100
0.150
0.200
0.250
0.300
5 10 15 20 25
....J| / cm3mol-1K
T / |J|
SPPA....
max
=18 and 19
SPPA....
max
=14 and 15
SPPA....
max
=10 and 11
SPPA....
max
=4 and 5
CSMC
SPPA....
max
=24 and 25
0.050
0.100
0.150
0.200
0.250
5 10 15 20 25
..|J| / cm3mol-1K
T / |J|
SPPA....
max
=18 and 19
SPPA....
max
=14 and 15
SPPA....
max
=10 and 11
SPPA....
max
=4 and 5
CSMC
(a) (b)
Fig. 6.26. χ|J | versus T/|J | plots as a function of the length of the spin correlation path for
two antiferromagnetic regular 2D networks: (a) square and (b) honeycomb. The results are
compared with the CSMC simulation (dots).
relation paths as a function of the correlation path length is shown. In Figure 6.26,
simulated χ|J | versus T/|J | curves are shown where the number of spin correlation
paths considered is limited by a prefixed maximum length of these paths (spin
path progressive addition method, SPPA).
These curves show better agreement with the MC simulation as this maximum
length increases, so in the infinite limit complete agreement is expected. However,
many of these spin correlation paths involve one or several loops. Thus, as has
been previously shown, these loops do not allow the factorisation of the partition
function and the total spin correlation function cannot be developed in individual
contributions, which invalidates this methodology [51]. Nevertheless, as these
closed paths become relevant only at low T/|J | values, then the simulated curves
reach the best agreement with the MC simulation at T/|J | > 8.7K. The threshold
T/|J | value for the applicability of the SPPA method is higher in a honeycomb
(6.6 K) than in a square networks, because the number of closed paths is lower
for a prefixed correlation path length in the first case. In Table 6.5, the temperature
expansion for Cur′ely’s law and for results obtained by SPPA, CSMC and

6.7 Some Complex Examples 217
Table 6.5. Coefficients of the temperature expansion series for 2D square and 2D honeycomb
networksobtainedbyCurely’slawandbySPPA,CSMCandHTEmethods[13,30,31,49,50].

Curely’s Law SPPA CSMC HTE
Square 2D
0 1.00000 1.00000 1.00000 1.00000
1 2.66667 2.66667 2.65737 2.66667
2 3.55555 5.33333 5.22918 5.33333
3 2.84444 9.95556 9.59006 9.95556
4 1.26420 17.69877 16.81410 16.90864
5 0.06020 31.24374 29.8237 27.24044
6 0.28896 53.99729 52.41480 42.21216
Honeycomb 2D
0 1.00000 1.00000 1.00000 1.00000
1 2.00000 2.00000 1.99966 2.00000
2 1.77778 2.66667 2.65477 2.66667
3 0.94815 3.02222 2.98714 3.02222
4 0.63210 3.31852 3.20519 3.31852
5 0.46655 3.67972 3.41378 3.67972
6 0.14448 3.83925 3.56714 3.57587

HTE methods are compared. In the CSMC method, the coefficients are obtained
from empirical laws presenting maximum terms β25 and β27 for the square and
honeycomb networks, respectively. These empirical laws are obtained from a fit
of the MC simulation data (see Section 6.3), which entails some uncertainties that
lead to very small discrepancies in the first coefficients of the expansion series
(see Table 6.5). A good agreement is obtained between the SPPA, CSMC and HTE
methods, whereas Cur′

ely’s law is certainly not efficient at describing the magnetic
behavior of 2D systems. The SPPA method diverges from the CSMC and the HTE
methods when the path length is long enough to consider loop diagrams. Differences found between CSMC and HTE methods are due to the limitations of this
last method at low T/|J|
values.

6.7 Some Complex Examples
The first example is a one-dimensional system with the formula
[{N(CH3)4}n][Mn2(N3)5(H2O)}n], which, from a magnetic point of view,
can be considered as a chain where there are magnetic couplings between near and
second neighbours (Figure 6.27) [53]. The interaction topology of this system has
been simplified by considering only two different exchange coupling constants. A


218 6 Monte Carlo Simulation: A Tool to Analyse Magnetic Properties

Fig. 6.27. Crystal structure and magnetic properties of [{N(CH3)4}n][Mn2(N3)5(H2O)}n]
[53]. The experimental data (circles) and the simulations by Monte Carlo methods (solid line)
are shown.

good agreement between the experimental and simulated data is obtained using
two different sets of parameters:
=2.001, J1 =1.57cm.1 and J2 =0.29cm.1
or
=2.007, J1 =0.66cm.1 and J2 =1.07cm.1. However, a careful analysis
of the structure (angles, bond lengths), as already done in Section 6.4, does not
reveal the set of parameters with the best physical meaning. Furthermore, the
determination of ferromagnetic interactions is always fairly inaccurate and the
use of a simplified interaction topology does not permit unambiguous assignment
[53].

Compound Csn[{Mn(N3)3}n] is a 3D solid where three different magnetic interactionsoccur(Figure6.28)[53].Regardingtheinteractiontopology,thissystemcan
be described as a stacking of alternating honeycomb planes. The magnetic behavior has been simulated using the set of parameters:
=2.029, J1 =0.76cm.1,
J2 =.4.3cm.1 and J3 =.3.3cm.1. The
values found for the interaction
througha μ-1,1-or μ-1,3-azidobridgearesimilartothosefoundinsimplersystems
[42, 53]. As in a previous example, the values of the J constants are corroborated by
the theoretical magneto-structural correlation performed by Ruiz et al. from DFT
calculations [44].

Fig. 6.28. Crystal structure and magnetic properties of Csn[{Mn(N3)3}n] [53]. The experimental data (circles) and the simulations by Monte Carlo methods (solid line) are shown.


6.7 Some Complex Examples 219
60
62
64
0 30 60
(a)
30
40
50
60
0 100 200 300
T / K
..
M
T / cm3 K mol-1
T / K
..
M
T
25
40
55
0 1400
(b)
T / K
..
M
T
J1
J2
J3
Fig. 6.29. Crystal structure, interaction topology and magnetic properties of
[Fe10Na2(O)6(OH)4(O2CPh)10(chp)6(H2O)2(MeCO)2] [54]. The experimental data
(circles) and the simulations by Monte Carlo methods (solid line) are shown.
Compound [Fe10Na2(O)6(OH)4(O2CPh)10(chp)6(H2O)2(MeCO)2] (chp = 6-
chloro-2-pyridonato) is an example of a high-nuclearity molecule (Figure 6.29)
[54]. This system is too big to be studied considering quantum spin moments, but
the CSMC method allows one to accurately simulate its magnetic behavior using
the coupling constant values: g = 2.0, J1 = .44 cm.1, J2 = .13 cm.1 and
J3 = .10 cm.1. Cano et al. rationalize the values of the the four coupling constants
taking into account the different bridging ligands, structural parameters and
some other data found in the literature. They conclude that the values found for the
four constants have a physical meaning.
Compound {[(tacn)6Fe8(μ3-O)2(μ2-OH)12]Br7(H2O)}Br·H2O (tacn = 1,4,7-
triazacyclononane) is one of a few examples of a single molecule magnet (Figure
6.30) [12, 55–57]. There are many interesting potential applications of these
systems. Although the study of the magnetic behavior of these systems is very
important, in some cases it is not yet possible to perform. This system is situated
at the limit where exact quantum solutions can be found. The simulated and experimental
χT versus T curves are shown in Figure 6.30. The theoretical curves
20.0
30.0
40.0
50.0
0 50 100 150 200 250 300
..T / cm3mol-1K
T / K
O
CSMC simulation
Quantum result
Experimetnal data
J1
J4
J3 J2
Fig. 6.30. Crystal structure, interaction topology and magnetic properties of {[(tacn)6Fe8(μ3-
O)2(μ2-OH)12]Br7(H2O)}Br·H2O [12, 55–57]. The experimental data (circles), the CSMC
simulation (solid line) and quantum solution (dashed line) are shown.

220 6 Monte Carlo Simulation: A Tool to Analyse Magnetic Properties

have been obtained from the exact quantum solutions and from the CSMC method
using, in both cases, the same parameter values. In the CSMC simulation, in order to better describe the magnetic behavior at very low temperatures, an extra
parameter (θ)
has been added to consider the magnetic intermolecular interaction
(
=2.0, J1 =.20cm.1, J2 =.120cm.1, J3 =.15cm.1, J4 =.35cm.1
and
=.2.2cm.1).

6.8 Conclusions and Future Prospects
In this chapter it has been established that the classical spin approach allows proper
analysisofthemagneticbehaviorofsystemswithhighlocalspinmoments (
≥2).
Nevertheless,thisapproachcannoteasilybeappliedtoagreatvarietyofsystems.In
these cases, it is possible to accomplish this objective using Monte Carlo methods,
which appears as a powerful tool in numerical integration to evaluate physical
properties. Thus, the Monte Carlo methods applied to a classical spin Heisenberg
model (CSMC) are able to study any system, whatever its complexity, and the
only limitation of this method is due to the classical spin approach. Unfortunately,
the simple CSMC method cannot be applied to systems that present small local
spin moments (
2). For such cases, it is possible to use alternative methods
although they are far more complicated. Among these methods are the Density
MatrixRenormalizationGroupandQuantumMonteCarlo.However,thisisanother
story, too long to be told in detail, and beyond the scope of the present chapter.