A Monte Carlo Wang-Landau_antiferromagnetic perovskite

[UFSJ]

Hi guys.

Using the Wang-Landau algorithm, I have tried to simulate the magnetic behavior of a bulk perovskite by the Monte Carlo method. I have used the JDoS description to build the Energy x magnetization phase space. By a referencing article, this magnetic system can be described by the Hamiltonian:
H = - 2J₁ ∑ S_i S_j - 2J₂ ∑ S_i
S_k - D ∑ S_i ², where J₁ is the exchange constant between two spins of the same basal plane (0.83 meV), J₂ is the exchange constant between two different basal planes (- 0.58 meV), D is the anisotropic constant (0.165 meV) and the experimental T_N = 139.5 K. This compound is an A-type Néel antiferromagnetic. Below, I share my Wang-Landau program to obtain the density of states (actually to obtain the ln[g(E,m)]) written in Fortran 90. Can anyone point out where my error is? I am obtaining a phase transition temperature at 75 K.

Thanks![CODE lang="fortran" title="Wang-Landau JDoS"]
PROGRAM Perov_4est_JDoS
IMPLICIT NONE

INTEGER, PARAMETER :: lz = 20 ! lattice length
INTEGER, PARAMETER :: ly = 20
INTEGER, PARAMETER :: lx = 20
INTEGER, PARAMETER :: n = lx*ly*lz ! number of spins
INTEGER, PARAMETER :: ff = 25 ! number of loops applying the modification factor e^0.5

REAL, PARAMETER :: flat = 0.95 ! the flatness parameter
INTEGER, PARAMETER :: cond = 1 ! if cond =1 the boundary conditions is ON, if cond = 0 is OFF

REAL, PARAMETER :: J1 = 2*0.83 ! exchange constant in the same basal plane in units of de meV
REAL, PARAMETER :: J2 = 2*(-0.58) ! exchange constant between basal planes in units of de meV
REAL, PARAMETER :: D = 0.1*J1 ! anisotropic constant in units of meV

REAL, PARAMETER :: E_max = n*((J1*16 - J2*8)/2 - D*4) ! maximum energy value

INTEGER, PARAMETER :: bin_E = 10**4 ! maximum number of different energy values to be get
INTEGER, PARAMETER :: bin_m = 10**4 ! maximum number of different magnetization values to be get

INTEGER :: r, k, q, r1, r2, k1, k2, px, py, pz, pV, pG, pr, pf, ping, pong, idum, sf
REAL :: Ei, Ef, mi, mf, E_latt, M_latt, E, m
DOUBLE PRECISION :: test, E_fund, h_sum, f, ln_Ohm_fund

INTEGER, DIMENSION (bin_E,bin_m) :: h, h_acum ! matrix Histograma and accumulated Histograma
DOUBLE PRECISION, DIMENSION (bin_E,bin_m) :: g ! matrix of the log of density of states
REAL, DIMENSION (bin_E) :: E_val ! array of energy values
REAL, DIMENSION (bin_m) :: m_val ! array of magnetization values

INTEGER, DIMENSION (n+1) :: s ! spin ! array os the values of the spin

INTEGER, DIMENSION (n) :: sv1 ! at right ! arrays of the neighbor's position along the spin chain
INTEGER, DIMENSION (n) :: sv2 ! at left
INTEGER, DIMENSION (n) :: sv3 ! in the next line
INTEGER, DIMENSION (n) :: sv4 ! in the bottom line
INTEGER, DIMENSION (n) :: sv5 ! in the up basal plane
INTEGER, DIMENSION (n) :: sv6 ! in the down basal plane

CHARACTER(LEN = 70) :: txt

REAL :: rd, xrand
idum = 123456789
s(n+1) = 0

! Array of first neighbors
DO pz=1,lz
DO py=1,ly
DO px=1,lx
r = px + lx*(py - 1) + lx*ly*(pz - 1)

sv1(r) = r + 1 ! at right
sv2(r) = r - 1 ! at left
sv3(r) = r + lx ! in the next line
sv4(r) = r - lx ! in the bottom line
sv5(r) = r + lx*ly ! up
sv6(r) = r - lx*ly ! down

END DO
END DO
END DO! Boundary condition
IF (cond==1) THEN ! boundary condition ON
DO pz=1,lz
DO py=1,ly
DO px=1,lx
r = px + lx*(py - 1) + lx*ly*(pz - 1)

IF(px==lx) THEN
sv1(r) = r - lx + 1 ! at right
END IF
IF(px==1) THEN
sv2(r) = r + lx - 1 ! at left
END IF

IF(py==ly) THEN
sv3(r) = r - lx*(ly - 1) ! in the next line
END IF
IF(py==1) THEN
sv4(r) = r + lx*(ly - 1) ! in the line before
END IFIF(pz==lz) THEN
sv5(r) = r - lx*ly*(lz - 1) ! up
END IF
IF(pz==1) THEN
sv6(r) = r + lx*ly*(lz - 1) ! down
END IF

END DO
END DO
END DOELSE IF (cond==0) THEN ! boundary condition OFF
DO pz=1,lz
DO py=1,ly
DO px=1,lx
r = px + lx*(py - 1) + lx*ly*(pz - 1)

IF(px==1) THEN
sv2(r) = n + 1 ! at left
END IF
IF(px==lx) THEN
sv1(r) = n + 1 ! at right
END IF

IF(py==1) THEN
sv4(r) = n + 1 ! in the line before
END IF
IF(py==ly) THEN
sv3(r) = n + 1 ! in the next line
END IF

IF(pz==1) THEN
sv5(r) = n + 1 ! up
END IF
IF(pz==lz) THEN
sv6(r) = n + 1 ! down
END IF

END DO
END DO
END DO

END IF

!!!!! Building the phase space E x m and updating the ln[g(E,m)] and histogram
DO r=1,bin_E
DO k=1,bin_m
g(r,k)= 1
h_acum(r,k) = 0
E_val(r) = 123456789 ! filling the energy and magnetization arrays with any wrong values. As a new state is observed, these values will be updating
m_val(k) = 123456789
END DO
END DO
DO k=1,n ! defining values for the spin-lattice
rd = xrand(idum)
IF (rd<=0.2) THEN
s(k) = - 2
ELSE IF ((0.2<rd).AND.(rd<=0.4)) THEN
s(k) = - 1
ELSE IF ((0.4<rd).AND.(rd<=0.6)) THEN
s(k) = 0
ELSE IF ((0.6<rd).AND.(rd<=0.8)) THEN
s(k) = 1
ELSE IF (0.8<rd) THEN
s(k) = 2
END IF
END DO

E_latt = 0
M_latt = 0
DO k=1,n ! obtaining the initial values of energy and magnetization of the chain
E_latt = E_latt - s(k)*( J1*( s(sv1(k)) + s(sv2(k)) + s(sv3(k)) + s(sv4(k)) )/2 + J2*( s(sv5(k)) + s(sv6(k)) )/2 + D*s(k) )
M_latt = M_latt + s(k)
END DO
Ei = E_latt
mi= M_latt
E_val(1) = Ei
m_val(1) = mi
DO pf=1,ff ! defining the value of the modification factor
f = exp(1.0)**( 1/(2**(pf - 1.0) ) )

DO r=1,bin_E
DO k=1,bin_m
h(r,k)= 0 ! histogram
END DO
END DO

print*, 'pf-', pf
ping = 0
DO WHILE (ping==0) ! the beginning of the loop responsible for the flatness test

DO pV=1,n*10**6 ! Monte Carlo loops testing new configuration of the chain

rd = xrand(idum) ! random choice of a spin "q" to be flipped
q = NINT(n*rd)
IF(q==0) THEN
q = n
END IF

rd = xrand(idum) ! Analyzing the flip of the spin "q"
IF (s(q)==-2) THEN
IF (rd<=0.25) THEN
sf = - 1
ELSE IF ((0.25<rd).AND.(rd<=0.5)) THEN
sf = 0
ELSE IF ((0.5<rd).AND.(rd<=0.75)) THEN
sf = 1
ELSE IF (0.75<rd) THEN
sf = 2
END IF

ELSE IF (s(q)==-1) THEN
IF (rd<=0.25) THEN
sf = - 2
ELSE IF ((0.25<rd).AND.(rd<=0.5)) THEN
sf = 0
ELSE IF ((0.5<rd).AND.(rd<=0.75)) THEN
sf = 1
ELSE IF (0.75<rd) THEN
sf = 2
END IF

ELSE IF (s(q)==0) THEN
IF (rd<=0.25) THEN
sf = - 2
ELSE IF ((0.25<rd).AND.(rd<=0.5)) THEN
sf = - 1
ELSE IF ((0.5<rd).AND.(rd<=0.75)) THEN
sf = 1
ELSE IF (0.75<rd) THEN
sf = 2
END IF

ELSE IF (s(q)==1) THEN
IF (rd<=0.25) THEN
sf = - 2
ELSE IF ((0.25<rd).AND.(rd<=0.5)) THEN
sf = - 1
ELSE IF ((0.5<rd).AND.(rd<=0.75)) THEN
sf = 0
ELSE IF (0.75<rd) THEN
sf = 2
END IF

ELSE IF (s(q)==2) THEN
IF (rd<=0.25) THEN
sf = - 2
ELSE IF ((0.25<rd).AND.(rd<=0.5)) THEN
sf = - 1
ELSE IF ((0.5<rd).AND.(rd<=0.75)) THEN
sf = 0
ELSE IF (0.75<rd) THEN
sf = 1
END IF

END IF
! Calculating the value of energy and magnetization of the possible final state with the flip of the spin "q"
Ef = Ei - (sf - s(q))*( J1*( s(sv1(q)) + s(sv2(q)) + s(sv3(q)) + s(sv4(q)) ) + J2*( s(sv5(q)) + s(sv6(q)) ) ) &
- D*(sf**2 - s(q)**2)
mf = mi + (sf - s(q))
DO r=1,bin_E ! identifying the position of the final energy along the Energy values array
IF (abs(E_val(r) - Ef) <= 0.001 ) THEN
r2 = r
GO TO 100
END IF
END DO

DO r=1,bin_E ! defining a position for the final energy along the Energy array, if it already not have been made
IF (abs(E_val(r) - 123456789) <= 0.001) THEN
E_val(r) = Ef
r2 = r
GO TO 100
END IF
END DO
100 CONTINUE

DO k=1,bin_m ! identifying the position of the final magnetization along the magnetization values array
IF (abs(m_val(k) - mf) <= 0.001 ) THEN
k2 = k
GO TO 200
END IF
END DO

DO k=1,bin_m ! defining a position for the final magnetization along the magnetization array, if it already not have been made
IF (abs(m_val(k) - 123456789) <= 0.001) THEN
m_val(k) = mf
k2 = k
GO TO 200
END IF
END DO
200 CONTINUE rd = xrand(idum) ! analyzing the possibility of changing from initial state to the final state and updating the histogram and density of states
IF (exp(g(r1,k1) - g(r2,k2)) >= 1) THEN
g(r2,k2) = g(r2,k2) + log(f)
h(r2,k2) = h(r2,k2) + 1
h_acum(r2,k2) = h_acum(r2,k2) + 1
s(q) = sf
r1 = r2
k1 = k2
ELSE IF (exp(g(r1,k1) - g(r2,k2)) >= rd) THEN
g(r2,k2) = g(r2,k2) + log(f)
h(r2,k2) = h(r2,k2) + 1
h_acum(r2,k2) = h_acum(r2,k2) + 1
s(q) = sf
r1 = r2
k1 = k2
ELSE ! Keeping in the initial state and updating the histogram and density of states
g(r1,k1) = g(r1,k1) + log(f)
h(r1,k1) = h(r1,k1) + 1
h_acum(r1,k1) = h_acum(r1,k1) + 1
END IFEND DO ! end of pV loop

!!!!!!!!!!!!!!! flatness test
h_sum = 0
pr = 0
DO r=1,bin_E
DO k=1,bin_m
IF (h_acum(r,k)>0) THEN ! considering just visited states, that is, the states with accumulated histogram not zero
pr = pr + 1
h_sum = h_sum + h(r,k)
END IF
END DO
END DO

pong = 0
DO r=1,bin_E
DO k=1,bin_m
IF (h_acum(r,k)>0) THEN
test = h(r,k)/(h_sum/pr)
IF (test < flat) THEN
pong = pong + 1
END IF
END IF
END DO
END DO
! if all considered states pass the flatness test, then "ping" changes to 1 and we go out of the flatness test, update the modification faction and repeat the process
IF (pong == 0) THEN
ping = 1
END IFEND DO ! final of the flatness loop

END DO ! final of the modification factor loop
!!!!!! Normalization of the density of states

ln_Ohm_fund = 0
E_fund = E_max ! starting considering the smallest value of energy as the maximum possible value
DO r=1,bin_E
DO k=1,bin_m
IF ((h_acum(r,k)>0).AND.(E_val(r) < E_fund)) THEN ! updating the smallest value of energy among the values computed
E_fund = E_val(r)
ln_Ohm_fund = g(r,k)
END IF
END DO
END DO

DO r=1,bin_E
DO k=1,bin_m
IF (h_acum(r,k)>0) THEN
g(r,k) = g(r,k) - ln_Ohm_fund + log(1.0) ! Updating the value of the density of states
END IF
END DO
END DO

! creating the output file with the energy and magnetization values and the respective ln[ g(E,m) ]
WRITE(txt,"(f4.2,f5.2,f4.2,I3,I3,I3,I0)") J1, J2, flat, ly, lx, lz, cond
WRITE(txt, "(A,f4.2,A,f5.2,A,f4.2,A,I3,A,I3,AI3,A,I0,A)") 'J1=',J1,'J2=',J2,'_ &
<h>=',flat,'lx=',lx,'ly=',ly,'lz=',lz,'cond=',cond,'.dat'
OPEN(2, FILE = trim(txt), STATUS = 'unknown')
DO r=1,bin_E
DO k=1,bin_m
IF (h_acum(r,k)>0) THEN
WRITE(2, *) E_val(r),';',m_val(k),';',g(r,k)
END IF
END DO
END DOEND PROGRAM

!cccccccccc Random function generator

FUNCTION XRAND(IDUM)
INTEGER IDUM
REAL XRAND
INTEGER, PARAMETER :: IM1=2147413563, IM2=2147413399, IMM1=IM1-1
INTEGER, PARAMETER :: IA1=40014,IA2=40692,IQ1=53661,IQ2=52774,IR1=12211,IR2=3791,NTAB=32,NDIV=1+IMM1/NTAB
REAL, PARAMETER :: AM=1./IM1,EPS=1.2E-7,RNMX=1.-EPS
INTEGER IDUM2,J,K,IV(NTAB),IY
SAVE IV,IY,IDUM2
DATA IDUM2/123456719/,IV/NTAB*0/,IY/0/
IF (IDUM .LE. 0)THEN
IDUM=MAX(-IDUM,1)
IDUM2=IDUM
DO J=NTAB+1,1,-1
K=IDUM/IQ1
IDUM=IA1*(IDUM-K*IQ1)-K*IR1
IF(IDUM .LT. 0)IDUM=IDUM+IM1
IF(J .LE. NTAB)IV(J)=IDUM
END DO
IY=IV(1)
END IF
K=IDUM/IQ1
IDUM=IA1*(IDUM-K*IQ1)-K*IR1
IF(IDUM .LT. 0)IDUM=IDUM+IM1
K=IDUM2/IQ2
IDUM2=IA2*(IDUM2-K*IQ2)-K*IR2
IF(IDUM2 .LT. 0)IDUM2=IDUM2+IM2
J= 1 + IY/NDIV
IY=IV(J)-IDUM2
IV(J)=IDUM
IF(IY .LT. 1)IY=IY+IMM1 XRAND=MIN(AM*IY,RNMX)
RETURN
END

[/CODE]

 

[berkeman]

Hi guys.

I have tried to simulate the magnetic behavior of a bulk perovskite by Monte Carlo method through the Wang-Landau algorithm. I have used the JDoS description of the phase space. By a reference article the magnetic system can be described by the hamiltonian:

Hi Professor, it looks like your post got cut off. Can you try again?

 

[UFSJ]

Hi Professor, it looks like your post got cut off. Can you try again?

Thank by the alert. And now, still cut off?

 

[berkeman]

Nope, looks fine now.