Monte Carlo Wang-Landau_antiferromagnetic perovskite

  • A
  • Thread starter UFSJ
  • Start date
  • #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

1. What is Monte Carlo Wang-Landau_antiferromagnetic perovskite?

Monte Carlo Wang-Landau_antiferromagnetic perovskite is a computational method used in statistical mechanics to study the properties of antiferromagnetic perovskite materials. It combines the Monte Carlo simulation technique with the Wang-Landau algorithm to efficiently sample the energy landscape of these materials.

2. How does Monte Carlo Wang-Landau_antiferromagnetic perovskite work?

The method works by randomly selecting a configuration of spins in the antiferromagnetic perovskite lattice and calculating its energy. The Wang-Landau algorithm then updates the energy histogram and adjusts the simulation parameters to ensure that all energy states are equally sampled. This process is repeated until the energy landscape is thoroughly explored and the density of states is accurately calculated.

3. What are the advantages of using Monte Carlo Wang-Landau_antiferromagnetic perovskite?

One of the main advantages is its ability to efficiently sample the energy landscape of complex antiferromagnetic perovskite materials, which can have a large number of possible spin configurations. It also allows for the calculation of thermodynamic properties, such as specific heat and magnetic susceptibility, which can provide valuable insights into the behavior of these materials.

4. What are some applications of Monte Carlo Wang-Landau_antiferromagnetic perovskite?

This method has been used to study a wide range of antiferromagnetic perovskite materials, including manganites, cuprates, and nickelates. It has also been applied to investigate phase transitions, magnetic ordering, and critical phenomena in these materials. Additionally, it has been used to design new materials with desired magnetic properties.

5. Are there any limitations to Monte Carlo Wang-Landau_antiferromagnetic perovskite?

One limitation is that the method requires a large number of simulation steps to accurately sample the energy landscape, which can be computationally expensive. It also assumes that the system is in thermal equilibrium, which may not always be the case for real materials. Additionally, it may not be suitable for studying materials with strong magnetic anisotropy or frustrated spin interactions.

Similar threads

Fortran Errors in Fortran code for solving Laplace's equation in 3D
    • Programming and Computer Science
Replies
4
Views
619
Fortran Need help with Jacobi relaxation method for Dirichlet boundary conditions
    • Programming and Computer Science
Replies
1
Views
945
I Prove series identity (Alternating reciprocal factorial sum)
    • General Math
Replies
4
Views
1K
Integration of Spherical Harmonics with a Gaussian (QM)
    • Calculus and Beyond Homework Help
Replies
1
Views
1K
A Matrices from Spherical Harmonics with Eigenvalue l+1
    • Quantum Physics
Replies
14
Views
1K
Testing floating point numbers for equality
    • Programming and Computer Science
Replies
7
Views
2K
Comp Sci Java: Monte Carlo experiment
    • Engineering and Comp Sci Homework Help
Replies
4
Views
2K
I Einstein brother-in-law elevator
    • Special and General Relativity
Replies
11
Views
196
I Deriving Lorentz Transformations Using Calculus
    • Special and General Relativity
Replies
18
Views
1K
[FDTD/FORTRAN] problem with tfsf boundary and berenger's pml
    • Atomic and Condensed Matter
Replies
8
Views
4K

Hot Threads

Recent Insights

Back
Top