A + *
The program is more general than the discussion as the
program implements two models,
A + * A2 + 2* The program is less general than the discussion as only a quadratic lattice is implemented.
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <assert.h>
#define NBRNUM 4
#define gascon 8.314
static struct SITE
{
struct SITE *nbr[NBRNUM];
int id;
int env;
int seq;
int flag;
} *site;
static struct QUEUE
{
double rate;
struct SITE *target;
int action;
} *queue;
static int count;
static double eqLimit;
static double omega;
static double desene;
static double difene;
static double preFacDif;
static double preFacDes;
static int size;
#define SQRT3 1.732050808
static double ran2(int offset)
{
unsigned int a = 16807;
unsigned int c = 0;
unsigned int m = 2147483647;
static int state = 1;
static int ready = 0;
if (offset > 0)
{
assert(ready == 0);
state = offset;
return 0.;
}
ready = 1;
state = (a * state + c) % m;
return ((double) state / m);
}
static void setEnvAtm(struct SITE * c)
{
int i;
c->env = 0;
for (i = 0; i < NBRNUM; i++)
if (c->nbr[i]->id > 0)
c->env++;
}
static void swapAtm(struct SITE * c, struct SITE * n)
{
int i, t;
if (c->id == 0)
{
if (n->id > 0)
{
for (i = 0; i < NBRNUM; i++)
c->nbr[i]->env++;
for (i = 0; i < NBRNUM; i++)
n->nbr[i]->env--;
}
}
else
{
if (n->id == 0)
{
for (i = 0; i < NBRNUM; i++)
c->nbr[i]->env--;
for (i = 0; i < NBRNUM; i++)
n->nbr[i]->env++;
}
}
t = c->id;
c->id = n->id;
n->id = t;
t = c->seq;
c->seq = n->seq;
n->seq = t;
}
static void swapNbr(struct SITE *c, int n)
{
swapAtm(c,c->nbr[n]);
}
static double calcene(void)
{
double e;
int i;
e = 0.;
for (i = 0; i < size * size; i++)
if (site[i].id > 0)
e += desene + omega * site[i].env;
return e;
}
static void walk(struct SITE * s, int difEq[], double *emin)
{
double e;
int i;
s->flag = s->env;
e = desene + 2. * omega * s->env;
if (e < *emin)
*emin = e;
if (difEq[s->env])
for (i = 0; i < NBRNUM; i++)
if (s->nbr[i]->id == 0 && s->nbr[i]->flag < 0)
{
swapAtm(s, s->nbr[i]);
walk(s->nbr[i], difEq, emin);
swapAtm(s, s->nbr[i]);
}
}
static void jump(struct SITE * s, int difEq[], double t)
{
double emin;
int i, i1, i2;
int last;
double level;
double sum;
double wghtSurf[NBRNUM + 1];
for (i = 0; i < size * size; i++)
site[i].flag = -1;
emin = desene + 2. * omega * s->env;
walk(s, difEq, &emin);
for (i = 0; i <= NBRNUM; i++)
if (desene + i * omega < emin - 5. * gascon * t)
wghtSurf[i] = 1.;
else if (desene + i * omega < emin + 5. * gascon * t)
wghtSurf[i] = exp(-(desene + 2. * i * omega - emin) / gascon / t);
else
wghtSurf[i] = 0.;
last = 0;sum=0.;
for (i = 0; i < size * size; i++)
if (site[i].flag >= 0)
{
queue[last].target = s;
queue[last].action = i;
assert(site[i].env >= 0 && site[i].env <= NBRNUM);
sum + =wghtSurf[site[i].flag];
queue[last].rate = sum;
last++;
}
if (last == 0)
return;
level = sum * ran2(0);
/* search by bisection */
i1 = -1;
i2 = last - 1;
while (i1+1 < i2)
{
i = (i1 + i2) / 2; assert(i>=0 && i<last);
if (level < queue[i].rate)
i2 = i;
else
i1 = i;
assert((i1 == -1 || queue[i1].rate <= level) && level < queue[i2].rate);
}
assert(i1+1 == i2 && (i2 == 0 || queue[i2 - 1].rate <= level) && level < queue[i2].rate);
swapAtm(s, &site[queue[i2].action]);
}
static void desOne(struct SITE * c)
{
int i;
c->id = 0;
count--;
for (i = 0; i < NBRNUM; i++)
setEnvAtm(c->nbr[i]);
}
static void desTwo(struct SITE * c, int n)
{
desOne(c);
desOne(c->nbr[n]);
}
static void output(int moves, int im)
{
int ia, is, ix, iy;
struct LIST
{
int id;
double x;
double y;
} *list;
if ((list = (struct LIST *) malloc((unsigned) size * size * sizeof(struct LIST))) == NULL)
{
(void) fprintf(stderr, "Not enough memory\n");
exit(1);
}
for (iy = 0; iy < size; iy++)
for (ix = 0; ix < size; ix++)
{
is = site[ix + size * iy].seq;
list[is].id = site[ix + size * iy].id;
list[is].x = (double) ix / (double) size;
list[is].y = (double) iy / (double) size;
}
if (im == 0)
{
(void) printf("\n");
(void) printf("2 %d 0 0 1\n", moves);
(void) printf("%d 0 0\n0 %d 0\n0 0 1\n", size, size);
(void) printf("%d %d\n", size * size);
for (ia = 0; ia < size * size; ia++)
{
(void) printf("%d", list[ia].id);
(void) printf(" X");
(void) printf(" %.3f %.3f %.3f", list[ia].x, list[ia].y, 0.);
(void) printf("\n");
}
}
else
{
(void) printf("%d 0 0\n0 %d 0\n0 0 1\n", size, size);
for (ia = 0; ia < size * size; ia++)
{
(void) printf("%d", list[ia].id);
(void) printf(" %.3f %.3f %.3f", list[ia].x, list[ia].y, 0.);
(void) printf("\n");
}
}
(void) free((char *) list);
}
static void equil(double temp,int moves)
{
double delene;
double oldene;
double newene;
int im, is, io;
oldene = calcene();
for (im = 0; im < moves; im++)
{
for (is = 0; is < size * size; is++)
if (site[is].id > 0)
{
io = NBRNUM * ran2(0);
if (site[is].nbr[io]->id == 0)
{
swapNbr(&site[is], io);
newene = calcene();
delene = newene - oldene;
if (delene < 0. || delene < 15. * gascon * temp && exp(-delene / gascon / temp) > ran2(0))
{
oldene = newene;
}
else
{
/* reject : swap id's back */
swapNbr(&site[is], io);
}
}
}
}
}
#define P(x,y) (x+y*size)
static void setup(void)
{
int j;
int i, ix, iy;
int xm, ym;
int xp, yp;
int t;
for (iy = 0; iy < size; iy++)
for (ix = 0; ix < size; ix++)
{
xm = (ix == 0) ? size - 1 : ix - 1;
ym = (iy == 0) ? size - 1 : iy - 1;
xp = (ix == size - 1) ? 0 : ix + 1;
yp = (iy == size - 1) ? 0 : iy + 1;
site[P(ix, iy)].seq = P(ix, iy);
site[P(ix, iy)].nbr[0] = &site[P(xp, iy)];
site[P(ix, iy)].nbr[1] = &site[P(ix, yp)];
site[P(ix, iy)].nbr[2] = &site[P(xm, iy)];
site[P(ix, iy)].nbr[3] = &site[P(ix, ym)];
}
for (i = 0; i < count; i++)
site[i].id = 1;
for (; i < size * size; i++)
site[i].id = 0;
for (i = 0; i < size * size; i++)
{
j = size * size * ran2(0);
if (site[i].id != site[j].id && ran2(0) > 0.5)
{
t = site[i].id;
site[i].id = site[j].id;
site[j].id = t;
}
}
for (i = 0; i < size * size; i++)
setEnvAtm(&site[i]);
}
static double rate1(double t)
{
int i, is;
double r;
double desRate[NBRNUM + 1];
r = 0.;
for (i = 0; i < NBRNUM; i++)
desRate[i] = preFacDes * exp(-(-desene - i * omega) / gascon / t);
for (is = 0; is < size * size; is++)
if (site[is].id > 0)
r += desRate[site[is].env];
return r;
}
static void model1(double tstart, double htrate, double tstep)
{
double difRate[NBRNUM + 1];
double desRate[NBRNUM + 1];
int difEq[NBRNUM + 1];
double dt;
int i, i1,i2,is;
int last;
double level;
double sum;
double t, tsum;
tsum = 0.;
for (;;)
{
t = tstart + htrate * tsum;
for (i = 0; i < NBRNUM; i++)
{
difRate[i] = preFacDif * exp(-(difene - i * omega) / gascon / t);
difEq[i] = (difRate[i] * tstep > eqLimit);
}
for (i = 0; i < NBRNUM; i++)
desRate[i] = preFacDes * exp(-(-desene - i * omega) / gascon / t);
for (is = 0; is < size * size; is++)
if (site[is].id > 0 && difEq[site[is].env])
jump(&site[is], difEq, t);
last = 0; sum=0.;
for (is = 0; is < size * size; is++)
if (site[is].id > 0)
{
if (!difEq[site[is].env])
for (i = 0; i < NBRNUM; i++)
if (site[is].nbr[i]->id == 0)
{
queue[last].action = i;
queue[last].target = &site[is];
assert(site[is].env >= 0 && site[is].env < NBRNUM);
sum += difRate[site[is].env];
queue[last].rate = sum;
last++;
}
queue[last].action = NBRNUM + 1;
queue[last].target = &site[is];
sum += desRate[site[is].env];
queue[last].rate = sum;
last++;
}
if (last == 0)
return;
dt = -log(ran2(0)) / sum;
if (tsum + dt * ran2(0) > tstep)
return;
level = sum * ran2(0);
/* search by bisection */
i1 = -1;
i2 = last - 1;
while (i1+1 < i2)
{
i = (i1 + i2) / 2; assert(i>=0 && i<last);
if (level < queue[i].rate)
i2 = i;
else
i1 = i;
assert((i1 == -1 || queue[i1].rate <= level) && level < queue[i2].rate);
}
assert(i1+1 == i2 && (i2 == 0 || queue[i2 - 1].rate <= level) && level < queue[i2].rate);
assert(queue[i2].target->id > 0);
if (queue[i2].action < NBRNUM)
swapNbr(queue[i2].target, queue[i2].action);
else
desOne(queue[i2].target);
tsum += dt;
}
}
static double rate2(double t)
{
int i, is;
double r;
double desRate[2 * NBRNUM + 1];
for (i = 0; i <= 2 * NBRNUM; i++)
desRate[i] = preFacDes * exp(-(2. * desene - i * omega) / gascon / t);
r = 0.;
for (is = 0; is < size * size; is++)
if (site[is].id > 0)
for (i = 0; i < NBRNUM / 2; i++)
if (site[is].nbr[i]->id == 1)
r += desRate[site[i].env + site[is].env];
return r;
}
static void model2(double tstart, double htrate, double tstep)
{
double difRate[NBRNUM + 1];
double desRate[2 * NBRNUM + 1];
int difEq[NBRNUM + 1];
double dt;
int i, i1,i2,is;
int last;
double rate;
double sum;
double t, tsum;
for (tsum = 0; tsum <= tstep; tsum += dt)
{
t = tstart + htrate * tsum;
for (i = 0; i <= NBRNUM; i++)
{
difRate[i] = preFacDif * exp(-(difene - i * omega) / gascon / t);
difEq[i] = (difRate[i] * tstep > eqLimit);
}
for (i = 0; i <= 2 * NBRNUM; i++)
desRate[i] = preFacDes * exp(-(2. * desene - i * omega) / gascon / t);
for (is = 0; is < size * size; is++)
if (site[is].id > 0 && difEq[site[is].env])
jump(&site[is], difEq, t);
last = 0;
sum=0.;
for (is = 0; is < size * size; is++)
if (site[is].id > 0)
{
if (!difEq[site[is].env])
for (i = 0; i < NBRNUM; i++)
if (site[is].nbr[i]->id == 0)
{
queue[last].action = i;
queue[last].target = &site[is];
assert(site[is].env >= 0 && site[is].env < NBRNUM);
sum += difRate[site[is].env];
queue[last].rate = sum;
last++;
}
for (i = 0; i < NBRNUM / 2; i++)
if (site[is].nbr[i]->id == 1)
{
queue[last].action = NBRNUM + i;
queue[last].target = &site[is];
assert(site[is].env >= 0 && site[is].env <= NBRNUM);
sum += desRate[site[i].env + site[is].env];
queue[last].rate = sum;
last++;
}
}
if (last == 0)
return;
dt = -log(ran2(0)) / sum;
if (tsum + dt * ran2(0) <= tstep)
{
rate = sum * ran2(0);
/* search by bisection */
i1 = -1;
i2 = last - 1;
while (i1+1 < i2)
{
i = (i1 + i2) / 2; assert(i>=0 && i<last);
if (rate < queue[i].rate)
i2 = i;
else
i1 = i;
assert((i1 == -1 || queue[i1].rate <= rate) && rate < queue[i2].rate);
}
assert(i1+1 == i2 && (i2 == 0 || queue[i2 - 1].rate <= rate) && rate < queue[i2].rate);
assert(queue[i2].target->id > 0);
if (queue[i2].action < NBRNUM)
swapNbr(queue[i2].target, queue[i2].action);
else
desTwo(queue[i2].target, queue[i2].action - NBRNUM);
}
}
}
int main(int argc, char *argv[])
{
int i;
int im;
int c;
void (*m) (double, double, double) = model1;
double (*r) (double) = rate1;
double tstep = 1.;
double htrate = 2.;
double conc = 0.95;
double tstart = 300.;
static int moves[2] = {0,500};
int seed = 0;
eqLimit = 10.;
omega = 0.;
preFacDif = 1.0e11;
preFacDes = 1.0e11;
desene = -60000.;
difene = 50000.;
size = 32;
while (--argc > 0 && (*++argv)[0] == '-')
{
c = *(++argv[0]);
switch (c)
{
case 'c':
if (--argc == 0
|| sscanf(*(++argv), "%lf", &conc) != 1
|| conc < 0. || conc > 1.)
{
(void) fprintf(stderr, "Could not read argument of option -%c\n", c);
return 1;
}
break;
case 'd':
if (--argc == 0
|| sscanf(*(++argv), "%lf", &preFacDif) != 1
|| preFacDif <= 0.
|| --argc == 0
|| sscanf(*(++argv), "%lf", &difene) != 1)
{
(void) fprintf(stderr, "Could not read argument of option -%c\n", c);
return 1;
}
break;
case 'e':
if (--argc == 0
|| sscanf(*(++argv), "%lf", &preFacDes) != 1
|| preFacDes <= 0.
|| --argc == 0
|| sscanf(*(++argv), "%lf", &desene) != 1)
{
(void) fprintf(stderr, "Could not read argument of option -%c\n", c);
return 1;
}
break;
case 'f':
if (--argc == 0
|| sscanf(*(++argv), "%lf", &eqLimit) != 1
||eqLimit < 1.)
{
(void) fprintf(stderr, "Could not read argument of option -%c\n", c);
return 1;
}
break;
case 'h':
if (--argc == 0
|| sscanf(*(++argv), "%lf", &tstep) != 1
|| tstep <= 0.)
{
(void) fprintf(stderr, "Could not read argument of option -%c\n", c);
return 1;
}
break;
case 'M':
if (--argc == 0
|| sscanf(*(++argv), "%d", &i) != 1
|| i < 1 || i > 2)
{
(void) fprintf(stderr, "Could not read argument of option -%c\n", c);
return 1;
}
if (i == 1)
{
m = model1;
r = rate1;
}
else
{
m = model2;
r = rate2;
}
break;
case 'm':
if (--argc == 0
|| sscanf(*(++argv), "%d", &moves[0]) != 1
|| moves[0] < 0
||--argc == 0
|| sscanf(*(++argv), "%d", &moves[1]) != 1
|| moves[1] < 1)
{
(void) fprintf(stderr, "Could not read argument of option -%c\n", c);
return 1;
}
break;
case 'o':
if (--argc == 0
|| sscanf(*(++argv), "%lf", &omega) != 1)
{
(void) fprintf(stderr, "Could not read argument of option -%c\n", c);
return 1;
}
break;
case 'r':
if (--argc == 0
|| sscanf(*(++argv), "%d", &seed) != 1
|| seed < 1)
{
(void) fprintf(stderr, "Could not read argument of option -%c\n", c);
return 1;
}
break;
case 's':
if (--argc == 0
|| sscanf(*(++argv), "%d", &size) != 1
|| size < 4)
{
(void) fprintf(stderr, "Could not read argument of option -%c\n", c);
return 1;
}
break;
case 't':
if (--argc == 0
|| sscanf(*(++argv), "%lf", &tstart) != 1
|| tstart <= 0.
|| --argc == 0
|| sscanf(*(++argv), "%lf", &htrate) != 1
|| htrate < 0.)
{
(void) fprintf(stderr, "Could not read argument of option -%c\n", c);
return 1;
}
break;
case '?':
(void) fprintf(stderr, "Usage: a.out [arguments]\n");
(void) fprintf(stderr, "\t-c <float> : start concentration\n");
(void) fprintf(stderr, "\t-d <float> <float> : preexponential factor and act energy for diffusion\n");
(void) fprintf(stderr, "\t-e <float> <float> : preexponential factor and act energy for desorption\n");
(void) fprintf(stderr, "\t-h <float> : timestep\n");
(void) fprintf(stderr, "\t-m <int> : number of timesteps\n");
(void) fprintf(stderr, "\t-M <int> : model\n");
(void) fprintf(stderr, "\t-o <int> : interaction between occupied sites\n");
(void) fprintf(stderr, "\t-r <int> : seed for RNG\n");
(void) fprintf(stderr, "\t-s <int> : size\n");
(void) fprintf(stderr, "\t-t <float> <float> : start temperature and heating rate\n");
default:
(void) fprintf(stderr, "-%c ?\n", c);
return 1;
}
}
if (seed > 0)
(void) ran2(seed);
count = conc * size * size;
if ((site = (struct SITE *) malloc((unsigned) size * size * sizeof(struct SITE))) == NULL
|| (queue = (struct QUEUE *) malloc((unsigned) (2 * NBRNUM) * count * sizeof(struct QUEUE))) == NULL)
{
(void) fprintf(stderr, "Not enough memory\n");
return 1;
}
setup();
for(im=0; im<moves[1]; im++)
{
equil(tstart,1);
(void) fprintf(stderr, "%d %.3f %.3f %.3f %.5f",
im, tstep * im,
tstart,
count / (double) (size * size),
0.);
if (count > 0)
(void) fprintf(stderr, " %.3f", calcene() / (count));
(void) fprintf(stderr, "\n");
}
output(moves[1],0);
for (im = 0; im < moves[1]; im++)
{
m(tstart + im * tstep * htrate, htrate, tstep);
(void) fprintf(stderr, "%d %.3f %.3f %.3f %.5f",
im, tstep * im,
tstart + htrate * tstep * im,
count / (double) (size * size),
r(tstart + (im + 1) * tstep * htrate) / (size * size));
if (count > 0)
(void) fprintf(stderr, " %.3f", calcene() / (count));
(void) fprintf(stderr, "\n");
output(moves[1],im + 1);
}
(void) free((char *) site);
(void) free((char *) queue);
return 0;
}