Monte Carlo Wang-Landau_antiferromagnetic perovskite

  • A
  • Thread starter UFSJ
  • Start date
In summary, the study on "Monte Carlo Wang-Landau antiferromagnetic perovskite" explores the application of the Wang-Landau algorithm in the context of antiferromagnetic materials, specifically perovskites. It focuses on the computational methods used to analyze the magnetic properties and phase transitions in these materials, highlighting the effectiveness of the Monte Carlo simulations in understanding complex magnetic behavior. The findings contribute to the knowledge of the thermodynamic characteristics of antiferromagnetic perovskites, which are significant for their potential applications in advanced electronic devices.
  • #1
UFSJ
15
2
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 = - 2J1 ∑ Si Sj - 2J2 ∑ Si
Sk - D ∑ Si 2, where J1 is the exchange constant between two spins of the same basal plane (0.83 meV), J2 is the exchange constant between two different basal planes (- 0.58 meV), D is the anisotropic constant (0.165 meV) and the experimental TN = 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
 
Last edited:
Physics news on Phys.org
  • #2
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?
 
  • #3
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?
 
  • Like
Likes berkeman
  • #4
Nope, looks fine now. :smile:
 
  • Like
Likes UFSJ

Similar threads

Fortran Errors in Fortran code for solving Laplace's equation in 3D
Replies
4
Views
2K
Fortran Need help with Jacobi relaxation method for Dirichlet boundary conditions
Replies
1
Views
1K
I Prove series identity (Alternating reciprocal factorial sum)
Replies
4
Views
2K
Integration of Spherical Harmonics with a Gaussian (QM)
Replies
1
Views
1K
Testing floating point numbers for equality
Replies
7
Views
2K
A Matrices from Spherical Harmonics with Eigenvalue l+1
Replies
14
Views
1K
[FDTD/FORTRAN] problem with tfsf boundary and berenger's pml
Replies
8
Views
4K
Comp Sci Why Is My Java Monte Carlo Simulation Not Approximating Pi Accurately?
Replies
4
Views
2K
I Einstein brother-in-law elevator
Replies
11
Views
967
I Deriving Lorentz Transformations Using Calculus
Replies
18
Views
2K

Hot Threads

Recent Insights

Back
Top