Program Heis1
*	molecular dynamics of Heisenberg chain using Runge-Kutta
*	i for lattice position, j for component
	include 'heis1.inc'
	real*8 save(3,5000),pe
	integer i,iorigin,itime,iconf 
	call initial
	do iconf = 1,ninitial
	   call newspin    ! new initial configuration
*	   equilibration
	   do i = 1,int(tequil/dt)
	      call runkut
	   enddo
	   do iorigin = 1,norigin
	      call torigin(save) ! saves the array in save
	      do itime = 0,tmax
	         call cfunct(save,itime) ! compute spin corr. funct.
	         do i = 1,int(0.1/dt)
	            call runkut
	         enddo
	      enddo
	   enddo
	enddo
	call printout
	stop
	end

	subroutine initial
	include 'heis1.inc'
	integer itime,ir
*	open command file and read the parameters
	open (unit = 1, file = 'heis1.com',status = 'old')
	read(1,*) dt
	read(1,*) n
	read(1,*) iseed
	read(1,*) tmax
	read(1,*) rmax
	read(1,*) norigin
	read(1,*) ninitial
* 	initialize correlation function
	do itime = 0,tmax
	   do ir = 0,rmax
	      C(itime,ir) = 0.0
	   enddo
	enddo
	end

	subroutine newspin
	include 'heis1.inc'
	integer i
	real*8 length
*	initial array of spins
	do i = 1,n
	   spin(1,i) = 2*rand(iseed)-1 
	   spin(2,i) = 2*rand(iseed)-1
	   spin(3,i) = 2*rand(iseed)-1   
	   length = (spin(1,i)**2+spin(2,i)**2+spin(3,i)**2)**0.5 
	   spin(1,i) = spin(1,i)/length
	   spin(2,i) = spin(2,i)/length
	   spin(3,i) = spin(3,i)/length
	enddo
	end

	subroutine torigin(save)
	include 'heis1.inc'
	real*8 save(3,5000) 
	integer i,j
	do i = 1,n
	   do j = 1,3
	      save(j,i) = spin(j,i)
	   enddo
	enddo
	end

	subroutine cfunct(save,itime)
	include 'heis1.inc'
	real*8 save(3,5000)
	integer i,j,ir,itime,icomp
	do ir = 0,rmax
	   do i = 1,n     ! location on saved configuration
	      do j = 1,3
*	         location of spin on actual lattice ir apart
	         icomp = (i+ir) - ((i+ir-1)/n)*n
	         C(itime,ir) = C(itime,ir) + spin(j,icomp)*save(j,i)
	      enddo
	   enddo
	enddo
	end

	subroutine runkut
	include 'heis1.inc'
	real*8 k,temp,product,constant
	integer i,j,m
	dimension k(3,5000,4),temp(3,5000),product(3)
	do i = 1,n
	   do j = 1,3
	      temp(j,i) = spin(j,i)
	   enddo
	enddo
	do m = 1,3
	   if(m.eq.3) then
	      constant = 1.0
	   else
	      constant = 0.5
	   endif
	   do i = 1,n
*	      compute first derivative for spin i
	      call deriv(i,temp,product)
	      do j = 1,3
	         k(j,i,m) = dt*product(j)
	      enddo
	   enddo
	   do i = 1,n
	      do j = 1,3
	         temp(j,i) = spin(j,i) + k(j,i,m)*constant
	      enddo
	   enddo
	enddo
	do i = 1,n
	   call deriv(i,temp,product)
	   do j = 1,3
	      k(j,i,4) = dt*product(j)
	   enddo
	enddo
	do i = 1,n
	   do j = 1,3
	      spin(j,i) = spin(j,i) + k(j,i,1)/6.0+k(j,i,2)/3.0
	      spin(j,i) = spin(j,i) + k(j,i,3)/3.0+k(j,i,4)/6.0
	   enddo
	enddo
	end

	subroutine deriv(i,temp,product)
*	compute S_i \times [S_{i+1} + S_(i-1}]
	include 'heis1.inc'
	real*8 sumv(3),product(3),leftv(3),rightv(3),spinv(3)
	real*8 temp(3,5000)
	integer i,j,left,right
*	store spin S_i in array spinv(j)
	do j = 1,3
	   spinv(j) = temp(j,i)
	enddo
	call pbc(i,left,right)
	do j = 1,3
	   leftv(j) = temp(j,left)
	   rightv(j) = temp(j,right)
	enddo
*	add vectors
	do j = 1,3
	   sumv(j) = leftv(j) + rightv(j)
	enddo
*	vector multiply vectors
	call vprod(spinv,sumv,product)
	end

	subroutine pbc(i,left,right)
	include 'heis1.inc'
	integer i,left,right
	left = i - 1
	right = i + 1
	if (i .eq. 1) then
	   left = n
	endif
	if (i .eq. n) then
	   right = 1
	endif
	end

	subroutine vprod(spinv,sumv,product)
	real*8 spinv(3),sumv(3),product(3)
*	three components of cross product of S_i with neighbors
	product(1) = spinv(2)*sumv(3) - spinv(3)*sumv(2)
	product(2) = spinv(3)*sumv(1) - spinv(1)*sumv(3)
	product(3) = spinv(1)*sumv(2) - spinv(2)*sumv(1)
	end

	subroutine energy(pe)
	include 'heis1.inc'
	real*8 pe
	integer i,j
	pe = 0.0
	do i = 1,n-1
	   do j = 1,3
	      pe = pe - spin(j,i)*spin(j,i+1)
	   enddo
	enddo
	do j = 1,3
	   pe = pe - spin(j,n)*spin(j,1)
	enddo
	end

	subroutine printout
	include 'heis1.inc'
*	print spin autocorrelation function
	integer weightfactor ! number of times C function evaluated
	weigthfactor = n*norigin*ninitial  ! number times correl fn calculated
	do ir = 0,rmax
	   write(*,*) ir
	   do itime = 0,tmax
	      write(*,*) C(itime,ir)/real(weightfactor)
	   enddo
	enddo
	end