chem/KPP/kpp/kpp-2.1/int/gillespie.f90

   1 MODULE KPP_ROOT_Integrator
   2
   3   USE KPP_ROOT_Parameters, ONLY : NVAR, NFIX, NREACT
   4   USE KPP_ROOT_Global, ONLY : TIME, RCONST, Volume
   5   USE KPP_ROOT_Stoichiom
   6   USE KPP_ROOT_Stochastic
   7   USE KPP_ROOT_Rates
   8   IMPLICIT NONE
   9
  10 CONTAINS
  11
  12 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  13   SUBROUTINE Gillespie(Nevents, T, SCT, NmlcV, NmlcF)
  14 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  15 !
  16 !   Gillespie stochastic integration
  17 !   INPUT:
  18 !       Nevents = no. of individual reaction events to be simulated
  19 !       SCT     = stochastic rate constants
  20 !       T       = time
  21 !       NmlcV, NmlcF = no. of molecules for variable and fixed species
  22 !   OUTPUT:
  23 !       T       = updated time (after Nevents reactions)
  24 !       NmlcV   = updated no. of molecules for variable species
  25 !
  26 !
  27 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  28       IMPLICIT NONE
  29       KPP_REAL:: T
  30       INTEGER :: Nevents
  31       INTEGER :: NmlcV(NVAR), NmlcF(NFIX)
  32       INTEGER :: i, m, issa
  33       REAL   :: r1, r2
  34       KPP_REAL :: A(NREACT), SCT(NREACT), x
  35
  36       DO issa = 1, Nevents
  37
  38           ! Uniformly distributed random numbers
  39           CALL RANDOM_NUMBER(r1)
  40           CALL RANDOM_NUMBER(r2)
  41           ! Avoid log of zero
  42           r2 = MAX(r2,1.e-14)
  43
  44           ! Propensity vector
  45           CALL  Propensity ( NmlcV, NmlcF, SCT, A )
  46           ! Cumulative sum of propensities
  47           DO i = 2, NREACT
  48             A(i) = A(i-1)+A(i);
  49           END DO
  50
  51           ! Index of next reaction
  52           x = r1*A(NREACT)
  53           DO i = 1, NREACT
  54             IF (A(i)>=x) THEN
  55               m = i;
  56               EXIT
  57             END IF
  58           END DO
  59
  60           ! Update time with time to next reaction
  61           T = T - LOG(r2)/A(NREACT);
  62
  63           ! Update state vector
  64           CALL MoleculeChange( m, NmlcV )
  65
  66      END DO
  67
  68   CONTAINS
  69
  70     SUBROUTINE PropensityTemplate( T, NmlcV, NmlcF, Prop )
  71       KPP_REAL, INTENT(IN)  :: T
  72       INTEGER, INTENT(IN)   :: NmlcV(NVAR), NmlcF(NFIX)
  73       KPP_REAL, INTENT(OUT) :: Prop(NREACT)
  74       KPP_REAL :: Tsave, SCT(NREACT)
  75     ! Update the stochastic reaction rates, which may be time dependent
  76       Tsave = TIME
  77       TIME = T
  78       CALL Update_RCONST()
  79       CALL StochasticRates( RCONST, Volume, SCT )
  80       CALL Propensity ( NmlcV, NmlcF, SCT, Prop )
  81       TIME = Tsave
  82     END SUBROUTINE PropensityTemplate
  83
  84   END SUBROUTINE Gillespie
  85
  86 END MODULE KPP_ROOT_Integrator