program Fortran_Figure7_6_b

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: totalmol=1000  ! number of molecules
 integer, parameter :: ranpom=1000000
 integer, parameter :: noboxes=40
 integer, parameter :: noiterations=100000000
 integer, parameter :: noburnin=10000000

 real(kind=long), parameter :: DiffConst=0.0001_long ! diffusion constant
 real(kind=long), parameter :: tmlom=0.02_long ! 2*maximallengtofmove
 real(kind=long), parameter :: repcon=0.0001_long ! 
 
 real(kind=long), dimension(totalmol) :: position
 real(kind=long), dimension(noboxes) :: distribution
 real(kind=long), dimension(ranpom) :: RanNumUnif
   
 integer :: in,iru,jj,qq,ss,ww
 real(kind=long) :: proposedposition,pomprob,potdiff
 real(kind=long) :: signalatpos,signalatproposedpos
 
 integer size, seed(2)
 data seed /6, 5/ 
      size = 2

 call random_seed(SIZE=size)
 call random_seed(PUT=seed(1:size))
 call random_number(RanNumUnif)
  
 iru=0
   
 do jj=1,totalmol
    iru=iru+1
    if (iru>ranpom) then
       call random_number(RanNumUnif)
       iru=1
    end if    
    position(jj)=RanNumUnif(iru)
 end do
  
 do jj=1,noboxes
    distribution(jj)=0
 end do
 
 !==============
 !burn in
 
 do in=1,noburnin

    ! choose randomly a molecule
       
    iru=iru+1
    if (iru>ranpom) then
       call random_number(RanNumUnif)
       iru=1
    end if    
    qq=int(real(totalmol)*RanNumUnif(iru))+1
    
    ! move randomly the chosen molecule
    
    iru=iru+1
    if (iru>ranpom) then
       call random_number(RanNumUnif)
       iru=1
    end if
    proposedposition=position(qq)+(RanNumUnif(iru)-0.5_long)*tmlom
    
    ! acceptance criterion
    
    if (position(qq)<0.5) then
       signalatpos=0.001_long*position(qq)
    else
       signalatpos=0.001_long-0.001_long*position(qq)
    end if      
    if (proposedposition<0.5) then
       signalatproposedpos=0.001_long*proposedposition
    else
       signalatproposedpos=0.001_long-0.001_long*proposedposition
    end if
    potdiff=-(signalatproposedpos-signalatpos)/DiffConst
    do ww=1,totalmol
       if (ww.NE.qq) then
          potdiff=potdiff-repcon/abs(position(qq)-position(ww))
          potdiff=potdiff+repcon/abs(proposedposition-position(ww))          
       end if
    end do
    if (potdiff>0.0_long) then
       pomprob=exp(-potdiff)
    else
       pomprob=1.0_long
    end if    
    if ((proposedposition<0.0).OR.(proposedposition>1.0)) then
       pomprob=0.0_long
    end if
    
    iru=iru+1
    if (iru>ranpom) then
       call random_number(RanNumUnif)
       iru=1
    end if
    if (RanNumUnif(iru)<pomprob) then
       position(qq)=proposedposition 
    end if
    
 end do
 
 ! ===============
 ! sampling
           
 do in=1,noiterations

    ! choose randomly a molecule
       
    iru=iru+1
    if (iru>ranpom) then
       call random_number(RanNumUnif)
       iru=1
    end if    
    qq=int(real(totalmol)*RanNumUnif(iru))+1
    
       
    ! move randomly the chosen molecule
    
    iru=iru+1
    if (iru>ranpom) then
       call random_number(RanNumUnif)
       iru=1
    end if
    proposedposition=position(qq)+(RanNumUnif(iru)-0.5_long)*tmlom
    
    ! acceptance criterion
    
    if (position(qq)<0.5) then
       signalatpos=0.001_long*position(qq)
    else
       signalatpos=0.001_long-0.001_long*position(qq)
    end if      
    if (proposedposition<0.5) then
       signalatproposedpos=0.001_long*proposedposition
    else
       signalatproposedpos=0.001_long-0.001_long*proposedposition
    end if   
    potdiff=-(signalatproposedpos-signalatpos)/DiffConst
    do ww=1,totalmol
       if (ww.NE.qq) then
          potdiff=potdiff-repcon/abs(position(qq)-position(ww))
          potdiff=potdiff+repcon/abs(proposedposition-position(ww))          
       end if
    end do     
    if (potdiff>0.0_long) then
       pomprob=exp(-potdiff)
    else
       pomprob=1.0_long
    end if    
    pomprob=exp(-potdiff)
    if ((proposedposition<0.0).OR.(proposedposition>1.0)) then
       pomprob=0.0_long
    end if
    
    iru=iru+1
    if (iru>ranpom) then
       call random_number(RanNumUnif)
       iru=1
    end if
    if (RanNumUnif(iru)<pomprob) then
       position(qq)=proposedposition 
    end if
    
    ! add to distribution
    
    do jj=1,totalmol
       ss=nint(real(noboxes)*position(jj)+0.5000_long)
       if (ss>noboxes) then
	   ss=noboxes
       end if	   
       if (ss<1) then
	   ss=1
       end if	   
       distribution(ss)=distribution(ss)+1.0_long
    end do   	         
 end do
 
 do ss=1,noboxes
    distribution(ss)=distribution(ss)/real(noiterations)
 end do   	         
     
 open(unit=1, file='data_Figure7_6_b.dat', status='replace')

100 format(f6.2)

       do ss=1,noboxes
          write(1,100) distribution(ss)
       end do   	    

 close(unit=1)

end program Fortran_Figure7_6_b