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

   1 void StochasticRates( double RCT[], double Volume, double SCT[] );
   2 void Propensity ( int V[], int F[], double SCT[], double A[] );
   3 void MoleculeChange ( int j, int NmlcV[] );
   4 double CellMass(double T);
   5 void Update_RCONST();
   6
   7 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
   8 void Gillespie(int Nevents, double Volume, double* T, int NmlcV[], int NmlcF[])
   9 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
  10 {
  11
  12       int i, m=0, event;
  13       double r1, r2;
  14       double A[NREACT], SCT[NREACT], x;
  15
  16      /* Compute the stochastic reaction rates */
  17       Update_RCONST();
  18       StochasticRates( RCONST, Volume, SCT );
  19
  20       for (event = 1; event <= Nevents; event++) {
  21
  22           /* Uniformly distributed ranfor (m numbers */
  23           r1 = (double)rand()/(double)RAND_MAX;
  24           r2 = (double)rand()/(double)RAND_MAX;
  25
  26           /* Avoid log of zero */
  27           r2 = (r2-1.0e-14) ? r2 : 1.0e-14;
  28
  29           /* Propensity vector */
  30           TIME = *T;
  31           Propensity ( NmlcV, NmlcF, SCT, A );
  32
  33           /* Cumulative sum of propensities */
  34           for (i=1; i<NREACT; i++)
  35             A[i] = A[i-1]+A[i];
  36
  37           /* Index of next reaction */
  38           x = r1*A[NREACT-1];
  39           for ( i = 0; i<NREACT; i++)
  40             if (A[i] >= x) {
  41               m = i+1;
  42               break;
  43             }
  44
  45           /* Update T with time to next reaction */
  46           *T = *T - log(r2)/A[NREACT-1];
  47
  48           /* Update state vector after reaction m */
  49           MoleculeChange( m, NmlcV );
  50
  51      } /* for event */
  52
  53 } /* Gillespie */