program Fortran_Figure7_6_a 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 :: nosavings=4 real(kind=long), parameter :: DiffConst=0.0001_long ! diffusion constant real(kind=long), parameter :: tmlom=0.02_long ! 2*maximallengtofmove real(kind=long), dimension(totalmol) :: position real(kind=long), dimension(noboxes,nosavings) :: distribution real(kind=long), dimension(ranpom) :: RanNumUnif integer :: in,iru,jj,qq,ss,noiterations,noburnin,ll real(kind=long) :: proposedposition,pomprob 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 noiterations=10000 noburnin=1000 do ll=1,nosavings noiterations=noiterations*10 noburnin=noburnin*10 write(*,*) ll,noiterations,noburnin do jj=1,noboxes distribution(jj,ll)=0.0_long end do 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 !============== !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 pomprob=exp((signalatproposedpos-signalatpos)/DiffConst) 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)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 pomprob=exp((signalatproposedpos-signalatpos)/DiffConst) 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)noboxes) then ss=noboxes end if if (ss<1) then ss=1 end if distribution(ss,ll)=distribution(ss,ll)+1.0_long end do end do do ss=1,noboxes distribution(ss,ll)=distribution(ss,ll)/real(noiterations) end do end do open(unit=1, file='data_Figure7_6_a.dat', status='replace') 100 format(f6.2) do ll=1,nosavings do ss=1,noboxes write(1,100) distribution(ss,ll) end do end do close(unit=1) end program Fortran_Figure7_6_a