[tsscds] Diff of /trunk/src/kmc.f90
Diff of /trunk/src/kmc.f90
 Parent Directory
|  Revision Log
|  Patch
| 8 |
! ndd species that dissapear |
! ndd species that dissapear |
| 9 |
character*80 :: title |
character*80 :: title |
| 10 |
integer, parameter :: dp = selected_real_kind(15,307) |
integer, parameter :: dp = selected_real_kind(15,307) |
| 11 |
integer :: m,nran,nr,nesp,inran,i,j,k,l,mu,ndisp,pd,kk,ijk,tcont,tint,lexit |
integer :: m,nran,nr,nesp,inran,i,j,k,l,mu,ndisp,pd,kk,ijk,tcont,tint |
| 12 |
real(dp) :: t,a0,r2a0,suma,ptot |
real(dp) :: t,a0,r2a0,suma,ptot |
| 13 |
real(dp) :: rnd(2) |
real(dp) :: rnd(2) |
| 14 |
integer,dimension(:),allocatable :: p,p0,re,pr,n,ndd,cont,pold,sig,sigold,rep |
integer,dimension(:),allocatable :: p,p0,re,pr,n,ndd,cont,pold,sig,sigold,rep |
| 40 |
print "(t3,a,20(i9))","Reacts:",re |
print "(t3,a,20(i9))","Reacts:",re |
| 41 |
print "(t3,a,20(i9))","Prods :",pr |
print "(t3,a,20(i9))","Prods :",pr |
| 42 |
print "(/,t3,a,1p,i10)","Step size =",tint |
print "(/,t3,a,1p,i10)","Step size =",tint |
| 43 |
big: do inran=1,nran |
outter: do inran=1,nran |
| 44 |
print "(/,t3,a,i4,/,t3,a,9999(i7))","Calculation number",inran,"Time(ps)",n |
print "(/,t3,a,i4,/,t3,a,9999(i8))","Calculation number",inran,"Time(ps)",n |
| 45 |
p=p0 |
p=p0 |
| 46 |
pold=p0 |
pold=p0 |
| 47 |
sig=0 |
sig=p0 |
| 48 |
t=0.d0 |
t=0.d0 |
| 49 |
tcont=0 |
tcont=0 |
|
lexit=0 |
|
| 50 |
rep=0 |
rep=0 |
| 51 |
print "(e10.4,9999(i7))",t,p |
print "(e10.4,9999(i8))",t,p |
| 52 |
do j=1,m |
do j=1,m |
| 53 |
a(j)=rate(j)*p(re(j)) |
a(j)=rate(j)*p(re(j)) |
| 54 |
enddo |
enddo |
| 55 |
a0=sum(a) |
a0=sum(a) |
| 56 |
do while(a0>0) |
inner: do while(a0>0) |
| 57 |
call random_number(rnd) |
call random_number(rnd) |
| 58 |
t=t-log(rnd(1))/a0 |
t=t-log(rnd(1))/a0 |
| 59 |
r2a0=rnd(2)*a0 |
r2a0=rnd(2)*a0 |
| 65 |
enddo s1 |
enddo s1 |
| 66 |
p(re(mu))=p(re(mu))-1 |
p(re(mu))=p(re(mu))-1 |
| 67 |
p(pr(mu))=p(pr(mu))+1 |
p(pr(mu))=p(pr(mu))+1 |
| 68 |
|
cont(mu)=cont(mu)+1 |
| 69 |
do j=1,m |
tcont=tcont+1 |
| 70 |
|
!Test of equilibrium |
| 71 |
|
eq: do j=1,m |
| 72 |
sigold(j)=sig(j) |
sigold(j)=sig(j) |
| 73 |
if(p(j)>pold(j)) sig(j)=0 |
if(p(j)>pold(j)) sig(j)=p(j) |
| 74 |
if(p(j)<pold(j)) sig(j)=1 |
if(p(j)<pold(j)) sig(j)=-p(j) |
| 75 |
if(sig(j).ne.sigold(j).and.p(j)>0.and.pold(j)>0) then |
if(sig(j)*sigold(j)<0) then |
| 76 |
rep(j)=rep(j)+1 |
rep(j)=rep(j)+1 |
| 77 |
if(rep(j)==1000) lexit=1 |
if(rep(j)==1000) exit inner |
| 78 |
endif |
endif |
| 79 |
enddo |
enddo eq |
| 80 |
if(lexit==1) exit |
!Test of equilibrium |
|
cont(mu)=cont(mu)+1 |
|
|
tcont=tcont+1 |
|
| 81 |
do j=1,m |
do j=1,m |
| 82 |
a(j)=rate(j)*p(re(j)) |
a(j)=rate(j)*p(re(j)) |
| 83 |
enddo |
enddo |
| 84 |
a0=sum(a) |
a0=sum(a) |
| 85 |
if(mod(tcont,tint)==0.or.a0==0) print "(e10.4,9999(i7))",t,p |
if(mod(tcont,tint)==0.or.a0==0) print "(e10.4,9999(i8))",t,p |
| 86 |
enddo |
enddo inner |
| 87 |
enddo big |
!Printing final populations |
| 88 |
print*,"Population of every species" |
print*,"Population of every species" |
| 89 |
do i=1,nesp |
do i=1,nesp |
| 90 |
print "(i6,i7)",i,p(i) |
print "(i6,i8)",i,p(i) |
| 91 |
enddo |
enddo |
| 92 |
|
!Printing statistics |
| 93 |
print*,"counts per process" |
print*,"counts per process" |
| 94 |
do i=1,m |
do i=1,m |
| 95 |
print "(i6,i20)",i,cont(i) |
print "(i6,i20)",i,cont(i) |
| 96 |
enddo |
enddo |
| 97 |
|
enddo outter |
| 98 |
end program kmc |
end program kmc |
|
|
Legend:
| Removed from v.289 |
|
| changed lines |
| |
Added in v.290 |
|
|
|
