Ising Model in One Dimension

c ising model  SI units employed ...linear chain
c The interaction described in the Ising model is only between nearest
c neighbor spins, and is -J for parallel spins and +J for antiparallel
c spins. The total energy can be expressed in the form
c E(ith particle)= -J0 * si* sj - mubohr*B*si,
c where si =  1. The first sum is over nearest neighbor
      real kb ,j0 ,mubohr
      dimension sold(0:50),snew(0:50),en(50)
      data T, B , kb,mubohr /1., 1. ,1.38e-23, 9.27e-24/
c nav ~no. of averages or iterations to be performed
      data nspin ,j0 ,iteran,nav/10, 1.72e-23 ,2  ,40/
c nseed integer n.e.1
      pi=2.*asin(1.)
      beta=1/(kb*T)
      sold(0)=0.
      sold(nspin+1)=0.
      snew(0)=0.
      snew(nspin+1)=0.
      z=.2
      do 20 i=1,nspin
      call random(pi,i,z,ran)
      z=ran
      sold(i)=-1.
c      if(i.le.5)sold(i)=1.
c      if(i.gt.5)sold(i)=-1.
c      if(ran.le.0.5)sold(i)=-1.
c      if(ran.gt.0.5)sold(i)=1.
      print*,'i,j,sold(i,j)=',i,sold(i)
20    continue
c first run
      do 30 i=1,nspin
      en(i)=-j0*(sold(i)*sold(i-1)+
     $ sold(i)*sold(i+1))  -mubohr*B*sold(i)
c      print*,'i,en(i)=',i,en(i)
30    continue
      eave=0.
      do 35 i=1,nspin
      eave=eave+en(i)
35    continue
      eave=eave/float(nspin)
      print*,'random initial eave=',eave
      print*,'   '
      do 300 k=1,nav
      do 40 i=1,nspin
      snew(i)=-sold(i)
      eflip=-j0*(snew(i)*sold(i-1)+
     $ snew(i)*sold(i+1))  -mubohr*B*snew(i)
      if(eflip.le.en(i))then
      en(i)=eflip
      sold(i)=snew(i)
      endif
      if(eflip.gt.en(i))then
      deltae=eflip-en(i)
      boltz=exp(-beta*deltae)
      call random(pi,iteran,z,ran)
      z=ran
      iteran=iteran+1
      if(ran.le.boltz)then
      en(i)=eflip
      sold(i)=snew(i)
      endif
      if(ran.gt.boltz)then
      endif
      endif
40    continue
      eave=0.
      do 135 i=1,nspin
      eave=eave+en(i)
135    continue
      eave=eave/float(nspin)
      print 100,k,eave
300   continue
100   format(2x,'iter, eave=',3x,i3,4x,e10.3)
      stop
      end

      subroutine random(pi,i,z,ran)
      x1=z*float(i)*pi
      y=sqrt(x1)
      z= y-int(y)
      ran=z
      return
      end





 i,j,sold(i,j)= 1 -1.
 i,j,sold(i,j)= 2 -1.
 i,j,sold(i,j)= 3 -1.
 i,j,sold(i,j)= 4 -1.
 i,j,sold(i,j)= 5 -1.
 i,j,sold(i,j)= 6 -1.
 i,j,sold(i,j)= 7 -1.
 i,j,sold(i,j)= 8 -1.
 i,j,sold(i,j)= 9 -1.
 i,j,sold(i,j)= 10 -1.
 random initial eave= -2.1690002E-023

  iter, eave=     1    -0.217E-22
  iter, eave=     2    -0.217E-22
  iter, eave=     3    -0.201E-22
  iter, eave=     4    -0.135E-22
  iter, eave=     5    -0.119E-22
  iter, eave=     6    -0.138E-22
  iter, eave=     7    -0.122E-22
  iter, eave=     8    -0.106E-22
  iter, eave=     9    -0.873E-23
  iter, eave=    10    -0.101E-22
  iter, eave=    11    -0.122E-22
  iter, eave=    12    -0.124E-22
  iter, eave=    13    -0.109E-22
  iter, eave=    14    -0.143E-22
  iter, eave=    15    -0.196E-22
  iter, eave=    16    -0.196E-22
  iter, eave=    17    -0.196E-22
  iter, eave=    18    -0.196E-22
  iter, eave=    19    -0.196E-22
  iter, eave=    20    -0.196E-22
  iter, eave=    21    -0.196E-22
  iter, eave=    22    -0.196E-22
  iter, eave=    23    -0.196E-22
  iter, eave=    24    -0.143E-22
  iter, eave=    25    -0.177E-22
  iter, eave=    26    -0.196E-22
  iter, eave=    27    -0.196E-22
  iter, eave=    28    -0.196E-22
  iter, eave=    29    -0.196E-22
  iter, eave=    30    -0.196E-22
  iter, eave=    31    -0.196E-22
  iter, eave=    32    -0.143E-22
  iter, eave=    33    -0.230E-22
  iter, eave=    34    -0.177E-22
  iter, eave=    35    -0.230E-22
  iter, eave=    36    -0.230E-22
  iter, eave=    37    -0.230E-22
  iter, eave=    38    -0.230E-22
  iter, eave=    39    -0.230E-22
  iter, eave=    40    -0.230E-22