Monte Carlo Wang-Landau_antiferromagnetic perovskite

UFSJ · Aug 28, 2023

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!

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

berkeman · Aug 28, 2023

UFSJ said:

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 · Aug 28, 2023

berkeman said:

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

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

berkeman · Aug 28, 2023

Nope, looks fine now.

Monte Carlo Wang-Landau_antiferromagnetic perovskite

Similar threads

Hot Threads

Recent Insights