View of /trunk/src/kmc.f90

program kmc
implicit none
! KMC computes the population vs time of a number of species involved in several dynamical processes
! rate: rate constant of a given processes  
! p:    population of a given species (p0 is its initial value)
! re (pr): reactant (product) for a given process
! ndisp number of species that dissapear
! ndd species that dissapear
character*80 :: title
integer, parameter :: dp = selected_real_kind(15,307)
integer :: m,nran,nr,nesp,inran,i,j,k,l,mu,ndisp,pd,kk,ijk,tcont,tint
real(dp) :: t,a0,r2a0,suma,ptot
real(dp) :: rnd(2)
integer,dimension(:),allocatable :: p,p0,re,pr,n,ndd,cont,pold,sig,sigold,rep
real (dp) ,dimension(:),allocatable :: a,rate
read(*,"(a80)") title
print "(t3,a80)",title
read(*,*) m,nesp,nran
allocate(re(m),pr(m),a(m),rate(m),p0(nesp),p(nesp),n(nesp),cont(m),pold(nesp),sig(nesp),sigold(nesp),rep(nesp))
n=(/ (l,l=1,nesp) /)
do i=1,m
read(*,*) rate(I),re(i),pr(i)
enddo
ptot=0.d0
print*, "at the beginning"
do ijk=1,m
   cont(ijk)=0
enddo
do i=1,nesp
read(*,*) p0(i)
print*, "P0(",i,")=",p0(i)
ptot=ptot+p0(i)
enddo
read(*,*) ndisp
allocate(ndd(ndisp))
if(ndisp>0) read(*,*) ndd
read(*,*) tint
print "(t3,a,i4,/,t3,a,i4,/,t3,a,i4,/)","# of calcs:",nran,"# of procs:",m,"# of specs:",nesp
print "(t3,a,1p,20(e9.2))","Rates:",rate
print "(t3,a,20(i9))","Reacts:",re
print "(t3,a,20(i9))","Prods :",pr
print "(/,t3,a,1p,i10)","Step size =",tint
outter: do inran=1,nran
   print "(/,t3,a,i4,/,t3,a,9999(i8))","Calculation number",inran,"Time(ps)",n
   p=p0
   pold=p0
   sig=p0
   t=0.d0
   tcont=0
   rep=0
   print "(e10.4,9999(i8))",t,p
   do j=1,m
      a(j)=rate(j)*p(re(j))
   enddo
   a0=sum(a)
   inner: do while(a0>0)
     call random_number(rnd)
     t=t-log(rnd(1))/a0      
     r2a0=rnd(2)*a0
     suma=0.d0
     pold=p
     s1: do mu=1,m
       suma=suma+a(mu)
       if(suma>=r2a0) exit s1
     enddo s1
     p(re(mu))=p(re(mu))-1
     p(pr(mu))=p(pr(mu))+1
     cont(mu)=cont(mu)+1
     tcont=tcont+1
!Test of equilibrium
     eq: do j=1,m
        sigold(j)=sig(j)
        if(p(j)>pold(j)) sig(j)=p(j)
        if(p(j)<pold(j)) sig(j)=-p(j)
        if(sig(j)*sigold(j)<0) then
           rep(j)=rep(j)+1
           if(rep(j)==1000) exit inner
        endif
     enddo eq
!Test of equilibrium
     do j=1,m
        a(j)=rate(j)*p(re(j))
     enddo
     a0=sum(a)
     if(mod(tcont,tint)==0.or.a0==0) print "(e10.4,9999(i8))",t,p 
   enddo inner
!Printing final populations
   print*,"Population of every species"
   do i=1,nesp
      print "(i6,i8)",i,p(i)
   enddo
!Printing statistics
   print*,"counts per process"
   do i=1,m
      print "(i6,i20)",i,cont(i)
   enddo
enddo outter
end program kmc