Finite temperatures


In this exercise we will investigate how statistics affect simulated measurables.

A single spin with an uniaxial anisotropy has a bi-stable, Ising like, magnetic state. At finite temperatures the stability of the magnetic state is not finite but follows an exponential (Arrhenius) relaxation behaviour. As seen in the lecture, ensemble averaging can be crucial for the analysis of such systems.

  • Investigate the amount of statistics that is needed to say something relevant about the life time of the magnetic state of the system.

  • Does the need of statistics change with system parameters? (temperature, anisotropy, external field)

  • Extra: Can you fit the relaxation behaviour to an Arrhenius function?

You can do a similar analysis for a finite 1d-chain by either modifying the single spin example, or starting from the SimpleSystems/HeisChain example.

  • Is there a difference by performing ensemble averaging compared to just increasing the system size?

  • Does the exchange interaction magnitude affect the stability of the spin chains?

An accessible article for those interested in spin chains and statistics can be found here: A. Vindigni Inorganica Chimica Acta, 361 3731 (2008).


In this exercise the thermalization rates in spin simulations will be investigated.

As mentioned in the lecture, thermalising the system before performing measurements is crucial for ensuring relevant results. Here we will investigate this for a simple cubic model system.

The initial inpsd.dat file looks as follows

simid simpcube
ncell     10        10       10                 System size            
BC        P         P         P                 Boundary conditions (0=vacuum, P=periodic)
cell      1.00000   0.00000   0.00000         
          0.00000   1.00000   0.00000
          0.00000   0.00000   1.00000
Sym       1                                     Symmetry of lattice (0 for no, 1 for cubic, 2 for 2d cubic, 3 for hexagonal)
do_prnstruct 1

posfile      ./posfile
momfile      ./momfile
exchange     ./jfile

Initmag     3                                   Initial config of moments (1=random, 2=cone, 3=spec., 4=file)

ip_mode     S
ip_mcanneal 1
            5000  300.0 1.00e-16 0.95

mode        S
Temp        300.00          K                    Temperature of the system
hfield      0.00000   0.00000   0.00000          Static H field

damping     0.01                                 Damping parameter (gamma)
nstep       10000                                Number of time-steps
timestep    1.000e-16       s                    The time step-size for the SDE-solver

plotenergy  1                                    Sample energies
do_avrg     Y                                    Measure averages
do_cumu     Y                                    Sample cumulative averages
do_spintemp Y                                    Measure spin temperature

and the almost trivial posfile and momfile are written as

1 1   0.000000  0.000000  0.000000
1 1 1.00 0.0 0.0 1.0 

The jfile, that will be changed during the exercise initially can look like

1 1  1.0       0.0       0.0       0.80000000000000
1 1  1.0       1.0       0.0       0.50000000000000

i.e. including nearest and next-nearest neighbours on the cubic lattice. Notice that since sym 1 is given in inpsd.dat, the jfile can be kept to a minimum of two lines.

  • Starting with the inputs as defined above, vary the simulation method and damping (where applicable) to investigate the thermalization rate of the system.

  • Is the thermalization faster when going from low to high temperatures or vice versa? Anything particular happening around Tc?

  • Change the sign of the next-nearest neighbour and redo the study. Is the magnetization a good measurable for determining the thermalization now?

Phase diagrams

Obtaining the M vs T relationship is probably the most common use case for Monte Carlo simulations on spin systems. In this exercise you can compare the MC functionalities of UppASD with a the ALPS package.

The system in question is here the 2d square lattice with NN exchange couplings. To compare with other model implementations this example uses the aunits Y flag which sets the temperature unit to the exchange strength J instead of Kelvin.

The inpsd.dat will here look as follows (to start with). Note the TEMP entries for initial and measurement temperatures.

simid  scHeis64 
ncell     64 64 1
BC        P         P         0                 Boundary conditions (0=vacuum,P=periodic)
cell      1.00000   0.00000   0.00000
          0.00000   1.00000   0.00000
          0.00000   0.00000   1.00000
Sym       1                                     Symmetry of lattice (0 for no, 1 for cubic, 2 for 2d cubic, 3 for hexagonal)

aunits Y
posfile    ./posfile
exchange   ./jfile
momfile    ./momfile
do_prnstruct 1

Mensemble 1
Initmag   3                                     (1=random, 2=cone, 3=spec., 4=file)

ip_mode   M                                     Initial phase parameters
ip_temp   TEMP                                  --
ip_mcNstep     5000                            --

mode      M                                     S=SD, M=MC
temp      TEMP                                  Measurement phase parameters
mcNstep   5000                                --

do_avrg   Y                                     Measure averages
plotenergy 1
do_cumu Y

and the posfile and momfile, and jfile files looks as

1   1   0.0 0.0 0.0                    
1 1 1.00000 0.0 0.0 1.0
 1 1    1.0    0.0    0.0    0.5000

Again, note that with aunits Y the exchange interaction in jfile is not defined in mRy but in the dimensionless energy scale of J (which is not Joule either).

In order to obtain the full M(T) curve, several simulations are needed at consecutive temperatures. This is preferably scripted, like in this example where we use a simple bash script


for Temp in 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0
    mkdir T$Temp/
    echo "Temp: " $Temp
    cp Base/* T$Temp/
    cd T$Temp/
    gsed -i "s/TEMP/$Temp/g" inpsd.dat
    ${SD_BINARY} > out.log
    cd ..

Here you either need to replace the ${SD_BINARY} expression, or export the location of your UppASD binary as the environment variable with the same name. You also need to use the same directory structure as intended, i.e. put the input files in a directory called Base and the script in the directory below.

  • Run the script and plot the resulting M(T) curve.

  • Compare with the reference data in the sc_64_ALPS.dat file

  • Are the simulation parameters “good enough” or are more thermalization/sampling steps needed to obtain an accurate M(T) curve?


Ground state determination is another situation where Monte Carlo and Atomistic spin dynamics simulations are used extensively. For isotropic ferromagnets this is usually a trivial task but for systems with exchange frustration and/or strongly anisotropic systems the ground state search can prove quite demanding.

In this exercise you will examine different kinds of magnetic systems with the task of finding their respective ground states.

Anisotropic systems

The first example starts from the square lattice of the previous exercise, but this time we will add a uniaxial anisotropy to the system. Copy the input files from the previous exercise and add a uniaxial anisotropy with a strength of 0.4J , where J still is the dimensionles unit of energy as earlier. An anisotropy is added by including the keyword anisotropy ./kfile to inpsd.dat and creating the anisotropy file kfile as per below. Also start with the jfile that only had positive exchange couplings.

1 1 -0.40   0.00   0.00000   0.00000   1.00000  0.0

Note that you can give all input files besides inpsd.dat more descriptive names if wanted.

  • With the anisotropy included, simulate the system starting from the ferromagnetic configuration by using initmag 3 and relax the system to a temperature approaching zero. You can visualize the structure using the GUI, or if it fails to run with the ASD_Tools/ script provided in the Github repository. <

This should give you a ferromagnetic order of the full sample, as expected.

  • Now change the starting guess to a fully disordered state by using initmag 1 and redo the relaxation. What state do you end up with?

The anisotropies in the system causes large energy barriers that can not be overcome at low temperatures. So unless the starting guess is good there is a risk that one ends up in a local minimum instead of the global one. To overcome these barriers one can employ higher temperatures which is typically done with the “Simulated annealing” (SA) approach, where one performs consecutive initial phase runs, starting from a high temperature and gradually reducing the temperature to the target value.

In UppASD, this is preferably done using the ip_mcanneal or ip_nphase constructs. An example is seen below.

ip_mode H
ip_mcanneal 5
1000 2.0
1000 0.5
1000 0.1
1000 0.01
1000 0.001

There are no global recipies for the optimal number if anneal sweeps and the temperature distibution so here empirical experience is your friend.

  • Try different SA sequences and see if you can relax the anisotropic system to its ferromagnetic ground state.

Dzyalosinskii-Moriya systems

Also without anisotropy, ground state relaxation can prove a formidable task.

simid SCsurf_T
ncell     128 128 1
BC        P         P         0                 Boundary conditions (0=vacuum,P=periodic)
cell      1.00000   0.00000   0.00000 
          0.00000   1.00000   0.00000
          0.00000   0.00000   1.00000

do_prnstruct 2
posfile   ./posfile
exchange  ./jfile
momfie    ./momfile
dm        ./dmfile

Initmag   1                                     (1=random, 2=cone, 3=spec., 4=file)

ip_mode N
ip_mcanneal 0

mode      S                                     M for MC and S for SD
temp      0.000
damping   0.500
Nstep     5000
timestep  1.000e-15       s                    The time step-size for the SDE-solver

plotenergy   1   
do_avrg      Y
do_cumu      Y

and the posfile, jfile, momfile as below

1 1 0.0 0.0 0.0
1 1  1.0       0.0       0.0      1.00000
1 1 -1.0       0.0       0.0      1.00000
1 1  0.0       1.0       0.0      1.00000
1 1  0.0      -1.0       0.0      1.00000
1 1 1.0 0.1 0.0 1.0 

Now we also have a Dzyaloshinskii-Moriya (DM) file dmfile that looks like

1 1  1.0000  0.0000 0.0000  0.81650  0.00000  0.00000
1 1 -1.0000  0.0000 0.0000 -0.81650 -0.00000 -0.00000
1 1  0.0000  1.0000 0.0000  0.00000  0.81650  0.00000
1 1  0.0000 -1.0000 0.0000 -0.00000 -0.81650 -0.00000

where the first five columns are the same as for a regular jfile i.e. atom sites and distance vector \(\mathbf{R}_{ij}\). The remaining three columns define the DM interaction vector \(\mathbf{D}_{ij}\). This system with ferromagnetic exchange and significant DM interactions has a spin-spiral ground state as shown below.


Fig. 1 Spin spiral ground state for a DM system.

  • For this setup, start from a random configuration with initmag 1 and relax the structure. What end state do you get?

  • If you do not get the spiral state, try to obtain it with simulated annealing. How close can you come? Is the energy stabilizing?

Glassy systems

As a final challenge, lets consider a spin glass system. A simple yet illustrative model for a spin glass system is given by the Edwards-Anderson model, where a nearest neighbour Hamiltonian on a cubic lattice, but with random exchange interactions are used.

The UppASD code can model Edwards-Anderson spin glasses by the keywords ea_model T and ea_sigma XX where ea_sigma contols the width of the Gaussian distribution of the randomized exchange interactions. Disregarding the ea_sigma keyword for now, we can define a simple cubic system with randomized exchange with the following inpsd.dat

simid bccFe100                                    
ncell     10        10       10                 System size            
BC        P         P         P                 Boundary conditions (0=vacuum, P=periodic)
cell      1.00000   0.00000   0.00000         
          0.00000   1.00000   0.00000
          0.00000   0.00000   1.00000
Sym       0
do_prnstruct 2

posfile    ./posfile
momfile    ./momfile
exchange   ./jfile

Initmag   1                                     Initial config of moments (1=random, 2=cone, 3=spec., 4=file)

ip_mode   N
ip_mcanneal 1
          10000 100.0 1.00e-16 0.95
          100 100.0 1.00e-16 0.95
          100 100.0 1.00e-16 0.95

mode      M
Temp      0.001         K                     Temperature of the system
hfield    0.00000   0.00000   0.00000           Static H field
mcnstep   10000                                Number of time-steps

do_avrg   Y                                     Measure averages
plotenergy 1
do_cumu Y

ea_model T

do_autocorr Y
acfile ./acfile

While the posfile and momfile are made as easy as possible

1 1   0.000000  0.000000  0.000000
1 1 1.00 1.0 0.0 0.0

Even though the exchange interactions will be randomized, we still have to set up a Hamiltonian with existing couplings since only the magnitude and not the directions will be randomized. Thus we use the following jfile

1 1  1.0       0.0       0.0       1.00
1 1 -1.0       0.0       0.0       1.00
1 1  0.0       1.0       0.0       1.00
1 1  0.0      -1.0       0.0       1.00
1 1  0.0       0.0       1.0       1.00
1 1  0.0       0.0      -1.0       1.00

With the glassy system setup. we end this exercise session with a friendly competition:

  • Without changing the system size and Hamiltonian, apply your thermalization and minimization skills to get the lowest possible energy for this system.

The participant with the lowest energy will get a symbolic price during the conference dinner.