Lecture Monte carlo simulations: Application to lattice models: Part III - TS. Ngô Văn Thanh

VSOP19, Quy Nhon 3-18/08/2013

Ngo Van Thanh, Institute of Physics, Hanoi, Vietnam.

Part III. Finite size effects and Reweighting methods

III.1. Finite size effects III.2. Single histogram method III.3. Multiple histogram method III.4. Wang-Landau method III.5. The applications

A Guide to Monte Carlo Simulations in Statistical Physics

Monte Carlo Simulation in Statistical Physics: An Introduction

D. Landau and K. Binder, (Cambridge University Press, 2009).

Understanding Molecular Simulation : From Algorithms to Applications

K. Binder and D. W. Heermann (Springer-Verlag Berlin Heidelberg, 2010).

Frustrated Spin Systems

D. Frenkel, (Academic Press, 2002).





Lecture notes PDF files : http://iop.vast.ac.vn/~nvthanh/cours/vsop/ Example code : http://iop.vast.ac.vn/~nvthanh/cours/vsop/code/

H. T. Diep, 2nd Ed. (World Scientific, 2013).

III.1. Finite size effects

 Consider a spin model on the lattice of finite size L

Using for determining the properties of the corresponding infinite system distinguish the order of the transition (first or second order transition) test the results of simulation of finite size systems Finite size scaling and critical exponents

Scaling forms of thermodynamic quantities



: are scaling functions



Scaling relations

(3.1)

(3.2)



The thermodynamic properties at the transition

reduce to proportionality constants



Scaling relations



(3.3)

If , then

The cumulants





fourth order cumulant of the order parameter



As the system size

(3.4)

For

For



fourth order cumulant of the energy



(3.5)

The logarithmic derivative of the nth power of the magnetization

For a continuous transition (second order transition)

(3.6)

where

(3.7)



Calculate the critical exponents

exponent :

exponent :

(3.8)

(3.9)

exponent :

(3.10)

exponent :

(3.11)

For a discontinuous transition (first order transition)

finite-size scaling laws



(3.12)

with



Ferromagnetic Ising spin model on simple cubic lattice

Ferromagnetic Ising spin model on simple cubic lattice

The behavior of VL for the q = 10 Potts model

in two dimensions

A Guide to Monte Carlo Simulations in Statistical Physics D. Landau and K. Binder

III.2. Single histogram method



Introduced by Ferrenberg and Swendsen (1988), Ferrenberg (1991) using histograms to extract information from Monte Carlo simulations Applying to the calculation of critical exponents

The theory Consider a Monte Carlo simulation performed at T = T0

  

generates system configurations with a frequency proportional to the Boltzmann weight, exp(-E/kBT). Histogram of energy and magnetization H(E, M) The probability of simultaneously observing the system

(3.13)

: is the number of configurations (density of states) with energy E and magnetization M

Partition function: (3.14)



over the range of

: provides an estimate for E and M values generated during the simulation

 We have

N is the number of measurements made

(3.15)



: is an estimate for the true density of states

(3.16)



From the histogram H(E, M) we can invert Eq. (3.15) to determine

the relationship between the histogram measured at K = K0 and the estimated probability distribution for arbitrary K

Replace in Eq. (3.13) with the expression for , and normalize the distribution.

(3.17)



The average value of any function of E and M :

Energy histogram vs energy for ferromagnetic Ising spin model on square lattice of size L = 64, at temperature K = 0.1540

(3.18)



Case study - one dimension single histogram From the histogram of energy H(E)



(3.19)



The estimated probability distribution for arbitrary K

The average value of any function of E



(3.20)

For example











The simulations Find the temperature T0 which is close to transition temperature Simulate the model by using the standard Monte Carlo simulation (Metropolis algorithm). Plot specific heat or susceptibility versus temperature, peaking the value of temperature at maxima of or . Choose the range of energy T0 = 2.26

Emin = -1.8, Emax = 1.0

 Do the simulation by using the single histogram method at T = T0

 - Gererating an array of N_BIN for energy range of interest

Program structure

- initialize the lattice

SUBROUTINE equilibrating ! Equilibrating process DO ieq = 1 to N_eq CALL monte_carlo_step() ENDDO END SUBROUTINE

! Averaging process SUBROUTINE averaging DO iav = 1 to N_av CALL monte_carlo_step() do-analysis ENDDO END SUBROUTINE

SUBROUTINE monte_carlo_step generate-one-sweep END SUBROUTINE

Source : single_histogram_ising.f90

PROGRAM SINGLE_HISTOGRAM_ISING IMPLICIT NONE INTEGER, PARAMETER :: n = 40,nms = n*n REAL,DIMENSION(n,n) :: s REAL,ALLOCATABLE,DIMENSION(:) :: e1,e2,m1,m2,he INTEGER :: i,j,ip,im,jp,jm,ie,ia,ms,neq,nav,ib,nbin REAL :: r,emin,emax,t,dt,temp,et1,et2 ! CALL RANDOM_SEED() neq = 200000 nav = 400000 t = 2.3; dt = 4.0 emax = -1.0*float(n*n); emin = -2.0*float(n*n) nbin = INT((emax-emin)/dt) + 1 ALLOCATE(he(nbin),e1(nbin),e2(nbin),m1(nbin),m2(nbin))

CALL spin_conf() ! initial configuration CALL equilibrating() he =0.0; e1 = 0.0; e2 = 0.0; m1 = 0.0; m2 = 0.0 CALL averaging() OPEN(UNIT=12,FILE='sh_ising.dat') i = 0 DO ib = 1,nbin IF (he(ib) > 0) THEN write(12,1) he(ib),e1(ib)/he(ib),e2(ib)/he(ib)& ,m1(ib)/he(ib),m2(ib)/he(ib) i = i + 1 ENDIF ENDDO 1 format(5(F16.6,1x)) CLOSE(12) OPEN(UNIT=12,FILE='sh_ising.info') write(12,*) t,i,n CLOSE(12) DEALLOCATE(he,e1,e2,m1,m2) ! CONTAINS ! list of subroutine/function

SUBROUTINE spin_conf ! CALL RANDOM_NUMBER(s) DO j = 1,n DO i = 1,n IF (s(i,j) > 0.5) THEN s(i,j) = 1.0 ELSE s(i,j) = -1.0 ENDIF ENDDO ENDDO END SUBROUTINE spin_conf ! SUBROUTINE equilibrating DO ie = 1, neq CALL monte_carlo_step() ENDDO END SUBROUTINE equilibrating

SUBROUTINE averaging DO ia = 1, nav CALL monte_carlo_step() temp = 0.0 DO i = 1,n ip = i + 1; im = i - 1 IF (i == n) ip = 1 IF (i == 1) im = n DO j = 1,n jp = j + 1; jm = j - 1 IF (j == n) jp = 1 IF (j == 1) jm = n temp = temp & - s(i,j)*(s(ip,j)+s(im,j)+s(i,jp)+s(i,jm)) ENDDO ENDDO temp = 0.5*temp ib = NINT((temp-emin)/dt) + 1

IF ((ib > 0) .AND. (ib <= nbin)) THEN he(ib) = he(ib) + 1.0 temp = temp/float(n*n) e1(ib) = e1(ib) + temp e2(ib) = e2(ib) + temp*temp ! temp = 0.0 DO j = 1,n DO i = 1,n temp = temp + S(i,j) ENDDO ENDDO temp = temp/float(n*n) m1(ib) = m1(ib)+ abs(temp) m2(ib) = m2(ib)+ temp*temp ENDIF ENDDO END SUBROUTINE averaging

SUBROUTINE monte_carlo_step DO ms = 1,nms CALL RANDOM_NUMBER(r) i = INT(r*float(n))+1 CALL RANDOM_NUMBER(r) j = INT(r*float(n))+1 ip = i + 1; im = i - 1 IF (i == n) ip = 1 IF (i == 1) im = n jp = j + 1; jm = j - 1 IF (j == n) jp = 1 IF (j == 1) jm = n et1 = -s(i,j)*(s(ip,j)+s(im,j)+s(i,jp)+s(i,jm)) et2 = -et1 CALL RANDOM_NUMBER(r) IF (r < exp(-(et2-et1)/t)) s(i,j) = -s(i,j) ENDDO END SUBROUTINE monte_carlo_step ! END PROGRAM SINGLE_HISTOGRAM_ISING

 Do analysis





Read the data from histogram MC simulation

Calculate the average value of physical quantities



Calculate the partition function

Source : single_histogram_analysis.f90

PROGRAM SINGLE_HISTOGRAM_ANALYSIS IMPLICIT NONE REAL,ALLOCATABLE,DIMENSION(:) :: e,e1,e2,m1,m2,he INTEGER :: n,ib,nbin,it,nt REAL :: r,tmin,tmax,t,t0,dt,z REAL :: et1,et2,mt1,mt2,cv,chi ! nt = 80 tmin = 2.1 tmax = 2.4 dt = (tmax-tmin)/float(nt-1) !

OPEN(UNIT=12,FILE='sh_ising.info') read(12,*) t0,nbin,n CLOSE(12) ALLOCATE(he(nbin),e(nbin),e1(nbin)& ,e2(nbin),m1(nbin),m2(nbin)) ! OPEN(UNIT=12,FILE='sh_ising.dat') DO ib = 1,nbin read(12,1) he(ib),e1(ib),e2(ib),m1(ib),m2(ib) ENDDO 1 format(5(F16.6,1x)) CLOSE(12) e = e1*float(n*n) ! OPEN(UNIT=12,FILE='ising_analysis.dat') t = tmin

DO it = 1,nt z = 0.0 DO ib = 1,nbin z = z + he(ib)*exp(e(ib)*(1/t0-1/t)) ENDDO et1 = 0.0; et2 = 0.0; mt1 = 0.0; mt2 = 0.0 DO ib = 1, nbin et1 = et1 + he(ib)*e1(ib)*exp(e(ib)*(1/t0-1/t)) et2 = et2 + he(ib)*e2(ib)*exp(e(ib)*(1/t0-1/t)) mt1 = mt1 + he(ib)*m1(ib)*exp(e(ib)*(1/t0-1/t)) mt2 = mt2 + he(ib)*m2(ib)*exp(e(ib)*(1/t0-1/t)) ENDDO et1 = et1/z; et2 = et2/z; mt1 = mt1/z; mt2 = mt2/z cv = float(n*n)*(et2-et1*et1)/t/t chi = float(n*n)*(mt2-mt1*mt1)/t write(12,1) t,et1,mt1,cv,chi t = t + dt ENDDO CLOSE(12) DEALLOCATE(he,e,e1,e2,m1,m2) END PROGRAM SINGLE_HISTOGRAM_ANALYSIS

 Ferromagnetic Ising spin model on the square lattice

Specific heat and susceptibility vs temperature with N = 40

red points : Standard MC result green points : Single histogram result

 Ferromagnetic Ising spin model on the square lattice

Specific heat and susceptibility vs temperature with N = 40, 50, 60

III.3. Multiple histogram method



Introduced by Ferrenberg and Swendsen (Phys. Rev. Lett. 63,1195, 1989) Developed from single histogram method The technique is known to reproduce with very high accuracy the critical exponents of second order phase transitions

The theory Consider n independent MC simulation





Each simulation at temperature Tj , the system has Nj configurations The overall probability distribution at temperature T

(3.21)

where

fi is chosen self-consistently using Eq. 3.21

(3.22)



The thermal average of a physical quantity A is then calculated by

Monte Carlo simulations









(3.23)



Perform n independent MC simulation



Choose the range of temperature [Tmin, Tmax] Divide the teperature into n points: Ti (i = 1, …, n) Choose the range of energy [Emin, Emax] Divide the energy into nbin



Perform the analysis







Read the simulation data Calculate self-consistently fi by iterating Eqs. (3.21) and (3.22) Do the calculations for the physical quantities by using (3.23)

Calculate the energy histogram, average energy, magnetization…

III.4. Wang-Landau method

Introduced by Wang and Landau, Phys. Rev. Lett. 86, 2050 (2001) It permits one to detect with efficiency weak first-order transitions Use to study classical statistical models with difficultly accessed microscopic states.

The theory

 Wang–Landau sampling



g(E) is defined as the number of spin configurations for any given E



The algorithm uses a random walk in energy space to obtain an accurate estimate for the density of states g(E)

(3.24)

g(E) independent of temperature, it can be used to find all properties of the system at all temperatures.

The classical partition function can either be written as a sum over all states or over all energies



 We begin with some simple ‘guess’ for the density of states: g(E) = 1 For improve g(E), spins are overturned according to the probability



(3.25)

Following each spin-flip trial the density of states is updateds

fi is a “modification factor” that is initially greater than 1 At the beginning of the random walk, the modification factor f can be as large as





(3.26)

A histogram H(E) records the number of times a state of energy E is visited Each time the energy histogram satisfies a certain ‘flatness’ criterion, f is reduced according to

and reset histogram H(E) = 0 for all E





The reduction process of the modification factor f is repeated several times until the final value ffinal is close enough to 1.0 The histogram is considered as flat if

x% is flatness criterion, it can be chosen between 70% and 95%

(3.27)



is the average (mean value) histogram.

The thermodynamic quantities can be evaluated by

(3.28)

(3.29)

where Z is the partition function defined by

(3.29)

The canonical distribution at a temperature T

(3.30)

Wang-Landau Monte Carlo scheme

(1) Set g(E) = 1; choose a modification factor (e.g. f0 = e1) (2) Choose an initial state (3) Choose a site i (4) Calculate the ratio of the density of states (5) Generate a random number r such that 0 < r < 1 (6) If r <  : then flip the spin

(7) Set

go to the next site and go to (4)

(8) If the histogram is not ‘flat’, (9) If the histogram is ‘flat’,

decrease f, e.g.

(10) Repeat steps (3)–(9) until (11) Calculate properties using final density of states g(E)

Source : wanglandau_ising.f90

Ferromagnetic Ising spin model on square lattice

PROGRAM WANG_LANDAU_ISING IMPLICIT NONE INTEGER, PARAMETER :: n = 10 REAL,DIMENSION(n,n) :: s REAL,ALLOCATABLE,DIMENSION(:) :: e1,e2,m1,m2,he,ge,pe INTEGER :: i,j,ip,im,jp,jm,ib,nbin,step,iold,inew REAL :: r,emin,emax,de,temp,eold,enew ! CALL RANDOM_SEED() emax = -0.5*float(n*n); emin = -2.0*float(n*n) + 16.0 de = 4.0 nbin = INT((emax-emin)/de) + 1 ALLOCATE(ge(nbin),he(nbin),e1(nbin)) ALLOCATE(pe(nbin),e2(nbin),m1(nbin),m2(nbin)) s = 1.0; s(1,1) = -1.0; s(5,5) = -1.0 e1(1) = emin DO ib = 2,nbin e1(ib) = e1(ib-1) + de ENDDO

pe =0.0; e2 = 0.0; m1 = 0.0; m2 = 0.0 CALL wanglandau() ! OPEN(UNIT=12,FILE='wl_ising.dat') temp = MINVAL(ge)-1.0 DO ib = 1,nbin IF (pe(ib) > 0.0) THEN e2(ib) = e2(ib)/pe(ib) m1(ib) = m1(ib)/pe(ib) m2(ib) = m2(ib)/pe(ib) ENDIF write(12,1) ge(ib)-temp,e1(ib),e2(ib),m1(ib),m2(ib) ENDDO 1 format(5(F16.6,1x)) CLOSE(12) OPEN(UNIT=12,FILE='wl_ising.info') write(12,*) nbin,n CLOSE(12) DEALLOCATE(ge,he,pe,e1,e2,m1,m2)

! list of functions/subroutines

CONTAINS SUBROUTINE wanglandau REAL :: logf,logffinal,xpc LOGICAL :: notfinish,notflat,notfilled xpc = 0.9; logf = 1.0; logffinal = 10.0**(-8.0) ge = 0.0; he = 0.0 eold = etot() iold = INT((eold-emin)/de) + 1 ge(iold) = logf; he(iold) = 1.0 step = 0 notfinish = .true. DO WHILE (notfinish) step = step + 1 notflat = .true. notfilled = .true. DO WHILE (notflat) CALL RANDOM_NUMBER(r) i = INT(r*float(n))+1 CALL RANDOM_NUMBER(r) j = INT(r*float(n))+1

! periodic boundary condition ip = i + 1; im = i - 1 IF (i == n) ip = 1 IF (i == 1) im = n jp = j + 1; jm = j - 1 IF (j == n) jp = 1 IF (j == 1) jm = n ! calculate the energy temp = -s(i,j)*(s(ip,j)+s(im,j)+s(i,jp)+s(i,jm)) enew = eold - 2.0*temp inew = INT((enew - emin)/de) + 1 IF (inew > 0 .AND. inew <=nbin) THEN CALL RANDOM_NUMBER(r) IF (ge(iold) - ge(inew) > log(r)) THEN s(i,j) = -s(i,j) iold = inew eold = enew ENDIF ENDIF

ge(iold) = ge(iold) + logf he(iold) = he(iold) + 1.0 IF (notfilled) THEN notfilled = ANY(he == 0.0) ELSE temp = xpc*sum(he)/float(nbin) IF (MINVAL(he) .GE. temp) notflat=.false. CALL averaging() ENDIF ENDDO print*, 'Step ',step,' is well done' logf = 0.5*logf he = 0.0 IF (logf < logffinal) notfinish = .false. ENDDO END SUBROUTINE wanglandau

REAL FUNCTION etot temp = 0.0 DO i = 1,n ip = i + 1; im = i - 1 IF (i == n) ip = 1 IF (i == 1) im = n DO j = 1,n jp = j + 1; jm = j - 1 IF (j == n) jp = 1 IF (j == 1) jm = n temp = temp & - s(i,j)*(s(ip,j)+s(im,j)+s(i,jp)+s(i,jm)) ENDDO ENDDO etot = 0.5*temp END FUNCTION etot

SUBROUTINE averaging ! pe(iold) = pe(iold) + 1.0 temp = eold/float(n*n) e2(iold) = e2(iold) + temp*temp temp = 0.0 DO j = 1,n DO i = 1,n temp = temp + S(i,j) ENDDO ENDDO temp = temp/float(n*n) m1(iold) = m1(iold)+ abs(temp) m2(iold) = m2(iold)+ temp*temp END SUBROUTINE averaging ! END PROGRAM WANG_LANDAU_ISING

Density of state vs energy for N = 10

 Wang Landau results for Ising case with N = 10

Energy and specific heat Magnetization and susceptibility

III.5. The applications Frustrated spin system

 A spin system is frustrated when one cannot find a configuration of spins to fully satisfy the interaction (bond) between every pair of spins In other words, the minimum of the total energy does not correspond to the minimum of each bond. This situation arises when:

the lattice geometry does not allow to satisfy all interaction bonds simultaneously there is a competition between different kinds of interactions acting on a spin by its neighbors

Effects of frustrated surface in Heisenberg thin films Study the effects of frustrated surfaces on the properties of thin films Lattice: made of stacked triangular layers Spin: Heisenberg spins with an Ising-like interaction anisotropy

 Methods:





 Model

 We consider a thin film made up by stacking Nz planes of triangular lattice of N × N lattice sites

V. T. Ngo and H. T. Diep, J. Appl. Phys. 91, 8399 (2002). V. T. Ngo and H. T. Diep, Phys. Rev. B 75, 035412 (2007).

Standard Monte Carlo simulations and Single histogram method Green’s function technique

 Hamiltonian







Interaction between two NN surface spins is equal to Js interaction between NN in interior layers are supposed to be ferromagnetic and all equal to J = 1 The two surfaces of the film are frustrated if Js is antiferromagnetic Js < 0

 MC simulation

The equilibrating time is about 106 MC steps per spin and the averaging time is 2 × 106 MC steps per spin. Films thickness Nz = 4 and plane size N = 24, 36, 48 and 60 Periodic boundary conditions are used in the XY planes

(3.25)

 Ground state





(diamonds) and (crosses) For

Angles between spins on layer 1 are all equal  Angles between vertical spins are 

 Monte Carlo results

Magnetization and susceptibility of first two layers for

Magnetization and susceptibility of first two layers for

Phase diagram in the space for

Phase I denotes the ordered phase with surface noncollinear spin configuration, Phase II indicates the collinear ordered state Phase III is the paramagnetic phase

and

 Layer susceptibilities versus T for L = 36, 48, 60 with Js = −0.5  Left (right) figure corresponds to the first- second-layer susceptibility



Critical exponents

 Maximum of surface-layer susceptibility versus L for L = 24, 36, 48, 60 with  The slope gives /

two-dimensional Ising model / = 1.75. For 3D, / = 1.97

Js = −0.5 (a,b) , Js = 0.5 (c) and I = | Is | = 0.1, in the ln-ln scale.

Critical behavior of magnetic thin films

 Methods:





Study the critical behavior of magnetic thin films as a function of the film thickness Lattice: simple cubic Spin: Ising

 Model

We consider a thin film made made from a ferromagnetic simple cubic lattice The size of the film is L × L × Nz





Standard Monte Carlo method Multiple histogram method

X.T. Pham Phu, V. T. Ngo and H. T. Diep, Surface Science 603 , 109–116 (2009).

Hamiltonian The periodic boundary conditions (PBC) is applied in the xy planes

 Monte Carlo results

Magnetization and susceptibility with



Susceptibility and V1 with



Two-dimensional Ising model  = 1

Critical exponents



Two-dimensional Ising model  = 1.75

Critical exponents



Critical exponent 3D = 0.613

Three-dimensional Ising model  = 0.63



Critical exponent 3D = 1.25

Three-dimension Ising model  = 1.24

Critical exponent  vs filmthickness

Critical exponent  vs filmthickness



From , we calculate the effective dimension of thin film

Stacked triangular antiferromagnets

 Methods:

Study the phase transition in the frustrated XY and Heisenberg spin model Lattice: Stacked triangular lattice (3D) Spin: XY and Heisenberg

  Wang-Landau method

 Model

We consider the stacking of triangular lattices in the z direction The system size is N × N × N Lattice size : N = 12, 18, 24, …, 120, 150



Standard Monte Carlo method

V. T. Ngo and H. T. Diep, J. Appl. Phys. 103, 07C712 (2008). V. T. Ngo and H. T. Diep, Phys. Rev. E 78, 031119 (2008).

Hamiltonian



XY antiferromagnets

Energy vs temperature for N = 84



XY antiferromagnets

magnetization vs temperature for N = 84



XY antiferromagnets

Energy histogram

 Heisenberg antiferromagnets

Energy and magnetization vs temperature for N = 120

 Heisenberg antiferromagnets

Energy histogram for N = 96 and 120

 Heisenberg antiferromagnets

Energy histogram for N = 150

Susceptibility vs temperature

 Heisenberg antiferromagnets

Maximum of susceptibility versus

Susceptibility vs temperature

N = 96, 108, 120, and 150 on a ln-ln scale

Crossover from first- to second-order transition in frustrated Ising antiferromagnetic films

 Methods:

Study the nature of this phase transition in the case of a thin film as a function of the film thickness Lattice: Face center cubic thin films Spin: Ising antiferromagnetic

  Wang-Landau method



Standard Monte Carlo method

 Model

Green’s function technique



Hamiltonian

X. T. Pham Phu, V. T. Ngo and H. T. Diep, Phys. Rev. E 79, 061106 (2009).

Consider a film of FCC lattice structure

 Ground state:



The spin configuration depends on the interface interaction Js

(a) ordering of type I for

(b) ordering of type II for

 The energy of a surface spin for two configurations



 Monte Carlo results for Js = J = -1 For the bulk case, L = Nz = 12

Energy vs temperature

Energy histogram with PBC (a) and without PBC (b) in z direction



For the Nz = 4

Energy histogram for L = 20, 30 and 40



Energy vs temperature

For the case Nz = 2 and L = 120

Energy histogram



Specific heat vs temperature

For the case Nz = 2 and L =120

Susceptibilities of sublattices (a) 1 and (b) 3

The latent heat E as a function of film thickness



Critical exponents for Nz = 2

The maximum value of V1

Maximum sublattice susceptibility

Fully frustrated simple cubic lattice

 Methods:

Study the nature of the phase transition in the fully frustrated Lattice: simple cubic lattice Spin: Ising, XY and Heisenberg

  Wang-Landau method

 Model



Standard Monte Carlo method

: for antiferromagnetic bond : for antiferromagnetic bond

V. T. Ngo, D. Tien Hoang and H. T. Diep, Phys. Rev. E 82, 041123 (2010) V. T. Ngo, D. Tien Hoang and H. T. Diep, Mod. Phys. Lett. 25, 929 (2011) V. T. Ngo, D. Tien Hoang and H. T. Diep, J. Phys.: Condens. Matter 23, 226002 (2011)

Hamiltonian

 Monte Carlo results for the case XY spin

Magnetization and susceptibility Energy and specific heat with N =24

 Monte Carlo results for the case XY spin

Energy histogram

The transition is clearly of first order at L = 36

Energy histogram for N =24

 Monte Carlo results for the case Heisenberg spin

Energy and specific heat Magnetization and susceptibility

 Monte Carlo results for the case Heisenberg spin

The transition is clearly of first order at L = 70

Energy histogram

 Monte Carlo results for the case Ising spin

Energy and specific heat

 Monte Carlo results for the case Ising spin

The transition is clearly of first order at L = 120 We conclude that a fully frustrated simple cubic lattice undergoes a first-order transition for Ising, XY and Heisenberg spin models

Energy histogram Maximum of the specific heat