A practical guide to kinetic Monte carlo simulations and classical molecular dynamics simulations by U. Burghaus, J. Stephan, L. vattuone, J.M. Rogowska NOVA Science Publisher, Inc., New York, 2006 ISBN 1-59454-531-6 Please note: 1) Include the header files (if they have not already been included). 2) Include int idum=(-1); // for the random number generator, include always this line in the examples if you would like to use the Numerical Recipes Random Number Generator. 3) Use the right rnd function (rnd1 or rnd4) in the program code, the one that is consistent with the library linked to the program. It does not make a difference which one you choose but do it consistently, please; or use one of the home made rnd functions, see the appendix. 4) If you take advantage of the "cut" and "past" commands of a text editor or word program then please note that text might be wraped around at the end of the text lines. This might cause problems if you try to compile the program. **************************************************************************************************** //////////////////////////////////////////////////////////////////////////////////////////////////// **************************************************************************************************** #include // standard "C" header #include // for simple mathematical functions #include // for simple input and output operations #include "nr.h" // this header defined the prototypes of functions and is not required for all compilers // for example void jacobn(float x, float y[], int n); // or float ran4(long *idum); #define NRANSI // some compiler require this statement, // which indicates that standard "C" code follows int idum=(-1); // for the random number generator, include always this line in the examples // use the right rnd function rnd1 or rnd4 - that does not make a difference // or use one of the home made rnd functions, see the appendix **************************************************************************************************** //////////////////////////////////////////////////////////////////////////////////////////////////// **************************************************************************************************** Listing 1: #include #include // standard "C" header #include // for simple mathematical functions #include // for simple input and output operations #include "nr.h" // from numerical recipes #define NRANSI // some compiler require this statement, // which indicates that standard "C" code follows int idum=(-1); // for the random number generator, include always this line in the examples //********** stochastic model of the //********** Langmuirian adsorption dynamics bgh // include the header lines etc. here (see ch. 2.2) //********** model parameters #define nx 50 // size of the lattice #define ny 50 int Average = 100; // average results of the last "Average" MC steps #define N 10000 // number of MC steps, i.e., total number of particles "exposed" on the surface float SaturationCoverage = 0.95; // iteration will be stopped if saturation is reached //********** program variables int Lattice[100][100]; // lattice float AdsorptionProbability[100]; // adsorption probability float Coverage[100]; // coverage void LangmuirianDynamics() { int i,j; // rows, columns index int Index; // results int x,y, x0,y0; // lattice site int Ndes,Nstick; // number of desorbing and adsorbing particles int Sum; // total number of adsorbed particles //********** initialize variables for (i=1; i<=nx+4; i++) for(j=1; j<=ny+4; j++) Lattice[i][j]=0; for(i=0;i0) Ndes++; // if lattice site filled then desorption else { Nstick++; // site free --> adsorption x0=x; y0=y ; // store adsorption site coordinates }; }; // averaging loop end if (Nstick==0) printf("increase averaging cycle number\n"); Lattice[x0][y0]=1; // adsorb the particle Sum++; // accumulates total number of adsorbed particles if ( ((int)(100*Sum) % (nx*ny)) == 0) { // show every 100th adsorbate AdsorptionProbability[Index]=(float)(Nstick)/(float)(Ndes+Nstick); Coverage[Index]=(float)(Sum)/(float)(nx*ny); printf("i=%d Sum=%d Nstick=%d Ndes=%d Coverage=%f S=%f\n", i,Sum,Nstick,Ndes,Coverage[Index],AdsorptionProbability[Index]); Index++; } }; // exposure loop end if ( (float)(Sum)/((float)(nx*ny)) < SaturationCoverage) printf("increase maximum number of exposed particles \n"); } int main(void) { LangmuirianDynamics(); } **************************************************************************************************** //////////////////////////////////////////////////////////////////////////////////////////////////// **************************************************************************************************** Listing 2: //********** stochastic model of the //********** Langmuirian adsorption dynamics bgh // include the header lines etc. here (see ch. 2.2) #include #include // standard "C" header #include // for simple mathematical functions #include // for simple input and output operations #include "nr.h" // from numerical recipes #define NRANSI // some compiler require this statement, // which indicates that standard "C" code follows int idum=(-1); // for the random number generator, include always this line in the examples //********** model parameters #define nx 50 // size of the lattice #define ny 50 int Average = 100; // average results of the last "Average" MC steps #define N 10000 // number of MC steps, i.e., total number of particles "exposed" on the surface float SaturationCoverage = 0.95; // iteration will be stopped if saturation is reached //********** program variables int Lattice[100][100]; // lattice float AdsorptionProbability[100]; // adsorption probability float Coverage[100]; // coverage void LangmuirianDynamics() { int i,j; // rows, columns index int Index; // results int x,y, x0,y0; // lattice site int Ndes,Nstick; // number of desorbing and adsorbing particles int Sum; // total number of adsorbed particles //********** initialize variables for (i=1; i<=nx+4; i++) for(j=1; j<=ny+4; j++) Lattice[i][j]=0; for(i=0;i0) Ndes++; // if lattice site filled then desorption else { Nstick++; // site free --> adsorption x0=x; y0=y ; // store adsorption site coordinates }; }; // averaging loop end if (Nstick==0) printf("increase averaging cycle number\n"); Lattice[x0][y0]=1; // adsorb the particle Sum++; // accumulates total number of adsorbed particles if ( ((int)(100*Sum) % (nx*ny)) == 0) { // show every 100th adsorbate AdsorptionProbability[Index]=(float)(Nstick)/(float)(Ndes+Nstick); Coverage[Index]=(float)(Sum)/(float)(nx*ny); printf("i=%d Sum=%d Nstick=%d Ndes=%d Coverage=%f S=%f\n", i,Sum,Nstick,Ndes,Coverage[Index],AdsorptionProbability[Index]); Index++; } }; // exposure loop end if ( (float)(Sum)/((float)(nx*ny)) < SaturationCoverage) printf("increase maximum number of exposed particles \n"); } int main(void) { LangmuirianDynamics(); } **************************************************************************************************** //////////////////////////////////////////////////////////////////////////////////////////////////// **************************************************************************************************** Listing 3-complete: //********** modified Kisliuk including defects and preadsorbates #include #include #include #include "nr.h" // from numerical recipes #define NRANSI //**********defect parameters : int def_label = 1; // = 2 pre-adsorbates // = 1 standard // = 0 include defects float def_cov = 0.1; // coverage of pre-exposed defects or preadsorbates float pdes_def = 0.5; // desorption probability from defect sites float pdes_def_prec = 0.5; // desorption probability for a thermalized particle from a precursor state float ed_def = 5.0; // heat of adsorption of defects/pre-adsorbates float st1_def = 0.0; // nn interaction energy between defects/pre-adsorbates and adsorbates float st1_def_def = 0.0; // nn interaction energies between pre-adsorbates and pre-adsorbates int f1_def; // number of 1st nearest neighbours int f2_def; // nu. of 2nd neighbour int f3_def; // nu. of 3rd neighbour int f1_def_ads; // number of 1st nearest neighbours int f2_def_ads; // nu. of 2nd neighbour int f3_def_ads; // nu. of 3rd neighbour //********** model parameters float pdes_clean = 0.9254; // probability of desorption from clean surface site // = 0.0 adsorbs always on clean site // = 1.0 desorbs always from clean site float pdes_occupied = 0.6; // probability of desorption from equilibrated state // = 1.0 desorbs always after equilibration // = 0.0 adsorbed always after equilibration float heat = -7.00; /* heat of adsorption */ float st1 = 0.0; /* 1st */ float nd2 = 0.0; /* 2nd */ float rd3 = 0.00; /* 3rd */ int max_hops = 150; /* max of extrinsic precursor hops */ int intrin_hop = 1; /* max of intrinsic precursor hops =1= no extrinsic precursor */ //********** controll parameters int nx = 40; /* size of the lattice */ int ny = 40; /* for tests choose 10x10 realistic results with 40x40 */ int grid_label=1; /* 1 = square lattice 0 = hexagonal latticse */ int thermalize_label=0; /* 0 = no thermalization of the configuration */ int thermalize_defects_label =0; /* 0 = no 1=yes */ int th_hops = 20; /* number of hops of each particle during thermalization */ int th_particle = 1600; /* number of particles which will be thermalized */ float stick_norm = -1.0; /* sticking probability will be normalized to ... if <0 no normalisation */ int average = 2000; /* average results of the last "average" monte carlo steps choose with respect to nx, ny 50x50 --> 20 10x10 --> 5 */ int cov_average = 5; float saturation_coverage = 0.98; /* iteration will be stopt if sat_cov...of CO is reached */ float energy_loop_preset = 300000; /* maximum energy loop number */ float stick_preset = 0.0; /* interation stops if S is lower then */ float des_loop_preset= 90000000; /* exposure loop must be larger than average especially in cases of Eij<0 */ //********** program variables int idum=(-1); /* for the random number generator */ float new_energy; /* energy of the new configuration */ float curr_energy; /* current energy of the configuration */ int ndes; /* number of desorption events per mc cycle, i.e. monte carlo steps necessary for a successfull adsorption event defined as one monte carlo cycle */ int nprec; /* number of hitting an already occupied site precursor steps */ int nstick; int boundary_flag=1; /* periodic boundary conditions y=1/n */ int sum; /* number of adsorbed particles */ int sum_def; /* number of preadsorbates/defects */ int index; /* storr data */ int equi_flag=0; /* equilibration yes/no */ float coverage_old, coverage_new; /* for the averaging of the coverages/S */ int energy_loop; /* counts the energy loops to prevent an invenite loop */ int lifetime; /* lifetime of the precursor state */ float averaged_lifetime; /* lateral interaction energies of adsorbates ***************/ int f1; /* number of 1st nearest neighbours */ int f2; /* nu. of 2nd neighbour */ int f3; /* nu. of 3rd neighbour */ float Stick[1000]; float coverage[1000]; int lattice [100][100]; /* adsorbate lattice */ int defects[100][100]; /* defect/preadsorbate lattice */ int temp_lattice[100][100]; /* for thermalized lattices */ float convig_energy[1000]; void periodic_boundary () { int i; for (i=3; i<(nx+2); i++) { lattice[1][i] = lattice[nx+1][i]; lattice[2][i] = lattice[nx+2][i]; lattice[nx+3][i] = lattice[3][i]; lattice[nx+4][i] = lattice[4][i]; lattice[i][1] = lattice[i][nx+1]; lattice[i][2] = lattice[i][nx+2]; lattice[i][nx+3] = lattice[i][3]; lattice[i][nx+4] = lattice[i][4]; }; } void calc_energy() // SQUARE lattice { int i, j; /* adsorbate - adsorbate interaction */ sum = f1 = f2 = f3 = 0; for (i=3; i<=nx+2; i++) { for (j=3; j<=ny+2; j++) { sum=sum+lattice[i][j]; f1 = ( lattice[i+1][j]+lattice[i][j+1]+ lattice[i-1][j]+lattice[i][j-1] ) * lattice[i][j] + f1; f2 = ( lattice[i+1][j+1]+lattice[i-1][j+1]+ lattice[i+1][j-1]+lattice[i-1][j-1] ) * lattice[i][j] + f2; f3 = ( lattice[i+2][j]+lattice[i][j+2]+ lattice[i-2][j]+lattice[i][j-2] ) * lattice[i][j] + f3; } } sum_def = f1_def = f2_def = f3_def = 0; f1_def_ads = f2_def_ads = f3_def_ads = 0; if (def_label==0) { /** preadsorbates - preadsorbates interaction*/ for (i=3; i<=nx+2; i++) { for (j=3; j<=ny+2; j++) { sum_def=sum_def+lattice[i][j]; f1_def = ( defects[i+1][j]+defects[i][j+1]+ defects[i-1][j]+defects[i][j-1] ) * defects[i][j] + f1_def; f2_def = ( defects[i+1][j+1]+defects[i-1][j+1]+ defects[i+1][j-1]+defects[i-1][j-1] ) * defects[i][j] + f2_def; f3_def = ( defects[i+2][j]+defects[i][j+2]+ defects[i-2][j]+defects[i][j-2] ) * defects[i][j] + f3_def; } } /** preadsorbates - adsorbates interaction*/ for (i=3; i<=nx+2; i++) { for (j=3; j<=ny+2; j++) { f1_def_ads = ( lattice[i+1][j]+lattice[i][j+1]+ lattice[i-1][j]+lattice[i][j-1] ) * defects[i][j] + f1_def_ads; f2_def_ads = ( lattice[i+1][j+1]+lattice[i-1][j+1]+ lattice[i+1][j-1]+lattice[i-1][j-1] ) * defects[i][j] + f2_def_ads; f3_def_ads = ( lattice[i+2][j]+lattice[i][j+2]+ lattice[i-2][j]+lattice[i][j-2] ) * defects[i][j] + f3_def_ads; } } }; /* if def_label */ new_energy=heat*sum+(st1*f1+nd2*f2+rd3*f3)/2 + /* adsorbates -- adsorbates interaction */ ed_def*sum_def + st1_def_def*f1_def/2.0 + /* preadsorbates -- preadsorbate interaction */ f1_def_ads*st1_def; /* preadsorbates -- adsorbates interaction */ } /*********************************************************************/ int calc_energy_hex() // HEXAGONAL lattice { int i, j; char answer[1]; /* adsorbate -- adsorbate interactions */ sum = f1 = f2 = f3 = 0; for (i=3; i<=nx+2; i++) { for (j=3; j<=ny+2; j++) { sum=sum+lattice[i][j]; f1 = ( lattice[i+1][j]+lattice[i-1][j]+ lattice[i][j+1]+lattice[i][j-1]+ lattice[i-1][j+1]+lattice[i+1][j-1]) * lattice[i][j] + f1; /* in that way both particles are counted as neighbours, so divide by 2 */ f2 = ( lattice[i-1][j+2]+lattice[i+1][j+1]+ lattice[i-1][j-1]+lattice[i+1][j-2]+ lattice[i+2][j-1]+lattice[i-2][j+1]) * lattice[i][j] + f2; f3 = ( lattice[i-2][j+2]+lattice[i][j+2]+ lattice[i+2][j]+lattice[i+2][j-2]+ lattice[i-2][j]+lattice[i-2][j]) * lattice[i][j] + f3; } } sum_def = f1_def = f2_def = f3_def = 0; f1_def_ads = f2_def_ads = f3_def_ads = 0; if (def_label==0) { /* defect -- defect interactions */ for (i=3; i<=nx+2; i++) { for (j=3; j<=ny+2; j++) { sum=sum+lattice[i][j]; f1_def = ( defects[i+1][j]+defects[i-1][j]+ defects[i][j+1]+defects[i][j-1]+ defects[i-1][j+1]+defects[i+1][j-1]) * defects[i][j] + f1_def; /* in that way both particles are counted as neighbours, so divide by 2 */ f2_def = ( defects[i-1][j+2]+defects[i+1][j+1]+ defects[i-1][j-1]+defects[i+1][j-2]+ defects[i+2][j-1]+defects[i-2][j+1]) * defects[i][j] + f2_def; f3_def = ( defects[i-2][j+2]+defects[i][j+2]+ defects[i+2][j]+defects[i+2][j-2]+ defects[i-2][j]+defects[i-2][j]) * defects[i][j] + f3_def; } } /* adsorbate -- defect interactions */ for (i=3; i<=nx+2; i++) { for (j=3; j<=ny+2; j++) { sum=sum+lattice[i][j]; f1_def_ads = ( lattice[i+1][j]+lattice[i-1][j]+ lattice[i][j+1]+lattice[i][j-1]+ lattice[i-1][j+1]+lattice[i+1][j-1]) * defects[i][j] + f1_def_ads; /* in that way both particles are counted as neighbours, so divide by 2 */ f2_def_ads = ( lattice[i-1][j+2]+lattice[i+1][j+1]+ lattice[i-1][j-1]+lattice[i+1][j-2]+ lattice[i+2][j-1]+lattice[i-2][j+1]) * defects[i][j] + f2_def_ads; f3_def_ads = ( lattice[i-2][j+2]+lattice[i][j+2]+ lattice[i+2][j]+lattice[i+2][j-2]+ lattice[i-2][j]+lattice[i-2][j]) * defects[i][j] + f3_def_ads; } } }; /* calculate total energy of the configuration */ new_energy=heat*sum+(st1*f1+nd2*f2+rd3*f3)/2 + /* adsorbates -- adsorbates interaction */ ed_def*sum_def + st1_def_def*f1_def/2.0 + /* preadsorbates -- preadsorbate interaction */ f1_def_ads*st1_def; /* preadsorbates -- adsorbates interaction */ return (0); } /*********************************************************************/ void check_boundary (int *x, int *y) /* coordinates oustide the lattice are reflected to the border of the lattice */ { if (*x<3) *x = nx+2; if (*x>(nx+2)) *x = 3; if (*y<3) *y = ny+2; if (*y>(ny+2)) *y = 3; } void thermalize () /* thermalize configuration by a random walk procedure */ { int i, j; /* old adsorption site */ int hop_direction; int x, y; /* new adsorption site */ char answer[10]; int k; /* number of hops for each particle*/ int l; /* number of prticle to thermalize */ /**********************************************/ /*****testing*****************************************/ /**********************************************/ int ff1; /* number of 1st nearest neighbours */ int ff2; /* nu. of 2nd neighbour */ int ff3; /* nu. of 3rd neighbour */ float test_sum; float E_before, E_after; if (coverage[index-1]>.05) { test_sum = ff1 = ff2 = ff3 = 0; for (i=3; i<=nx+2; i++) { for (j=3; j<=ny+2; j++) { sum=sum+lattice[i][j]; ff1 = ( lattice[i+1][j]+lattice[i][j+1]+ lattice[i-1][j]+lattice[i][j-1] ) * lattice[i][j] + ff1; /* in that way both particles are counted as neighbours, so divide by 2 */ ff2 = ( lattice[i+1][j+1]+lattice[i-1][j+1]+ lattice[i+1][j-1]+lattice[i-1][j-1] ) * lattice[i][j] + ff2; ff3 = ( lattice[i+2][j]+lattice[i][j+2]+ lattice[i-2][j]+lattice[i][j-2] ) * lattice[i][j] + ff3; } } E_before=heat*sum+(st1*ff1+nd2*ff2+rd3*ff3)/2; } /**********************************************/ /*****thermalizing*****************************/ /**********************************************/ //printf("before\n"); show_config(); for (l=1; l<=th_particle;l++) { /* choose "th_particle" to thermalize */ i = (int)((ran4(&idum)*nx)+3); /* choose particle to thermalize */ j = (int)((ran4(&idum)*ny)+3); if (lattice[i][j]==1) { /* hop only filled sites */ for (k=1; k=3) & (x<=(nx+2)) & (y>=3) & (y<=(ny+2)) & (lattice[x][y]==0)) { /* hop successfull */ lattice[x][y]=1; lattice[i][j]=0; i=x; j=y; }; }; /* for k */ }; /* if lattice */ }; /* for l */ //printf("after\n"); show_config(); /**********************************************/ /****testing******************************************/ /**********************************************/ if (coverage[index-1]>0.05) { test_sum = ff1 = ff2 = ff3 = 0; for (i=3; i<=nx+2; i++) { for (j=3; j<=ny+2; j++) { sum=sum+lattice[i][j]; ff1 = ( lattice[i+1][j]+lattice[i][j+1]+ lattice[i-1][j]+lattice[i][j-1] ) * lattice[i][j] + ff1; /* in that way both particles are counted as neighbours, so divide by 2 */ ff2 = ( lattice[i+1][j+1]+lattice[i-1][j+1]+ lattice[i+1][j-1]+lattice[i-1][j-1] ) * lattice[i][j] + ff2; ff3 = ( lattice[i+2][j]+lattice[i][j+2]+ lattice[i-2][j]+lattice[i][j-2] ) * lattice[i][j] + ff3; } } E_after=heat*sum+(st1*ff1+nd2*ff2+rd3*ff3)/2; if (E_after>=E_before) { printf("thermalizing error \n"); printf("E_after = %f E_before = %f \n",E_after,E_before); scanf("%s",answer); } } /**********************************************/ /**********************************************/ /**********************************************/ } /**********************************************************************/ void thermalize_defects () /* thermalize configuration by a random walk procedure */ { int i, j; /* old adsorption site */ int hop_direction; int x, y; /* new adsorption site */ char answer[10]; int k; /* number of hops for each particle*/ int l; /* number of prticle to thermalize */ for (l=1; l<=th_particle;l++) { /* choose "th_particle" to thermalize */ i = (int)((ran4(&idum)*nx)+3); /* choose particle to thermalize */ j = (int)((ran4(&idum)*ny)+3); if (defects[i][j]==1) { /* hop only filled sites */ for (k=1; k=3) & (x<=(nx+2)) & (y>=3) & (y<=(ny+2)) & (defects[x][y]==0)) { /* hop successfull */ defects[x][y]=1; defects[i][j]=0; i=x; j=y; }; }; /* for k */ }; /* if lattice */ }; /* for l */ } /*********************************************************************/ int make_a_hop(int hop_direction, int *xx, int *yy) /* x,y: position of the chosen neighbour site */ /* hop_direction: direction of the hop, expected to be chosen by random */ { char answer[10]; int x,y; /* printf("direction %d \n",hop_direction); */ x=*xx; y=*yy; /* use case statement */ /******************** 1st neighbor site hops ***********************/ if (hop_direction==1) { /* 1st top */ x=x-1; y=y; }; if (hop_direction==2) { /* 1st bottom */ x=x+1; y=y; }; if (hop_direction==3) { /* 1st right */ x=x; y=y+1; }; if (hop_direction==4) { /* 1st left */ x=x; y=y-1; }; /******* 2nd neigbhour hops **************************************/ if (hop_direction==5) { /* 2nd right top */ x=x-1; y=y+1; }; if (hop_direction==6) { /* 2nd right bottom */ x=x+1; y=y+1; }; if (hop_direction==7) { /* 2nd left top */ x=x-1; y=y-1; }; if (hop_direction==8){ /* 2nd left bottom */ x=x+1; y=y-1; }; check_boundary (&x, &y); *xx=x; *yy=y; return (0); } /********************************************************************/ void pre_expose() { int def_sum; int x,y; def_sum=0; while ( (def_sum<=(def_cov*nx*ny)) ) { x = (int)((ran4(&idum)*nx)+3); y = (int)((ran4(&idum)*ny)+3); if (defects[x][y]==0) { defects[x][y]=1; def_sum++; }; /* if defects */ }; /* while */ } void defects_mode(int *xx, int *yy, int *flag) { int x,y,x0,y0; int stick_flag; int extr_flag; float prob; int hop_direction; int hop_numbers; int l; x=*xx; y=*yy; stick_flag=*flag; //********** direct adsorption if ( (def_label==0) & (lattice[x][y]==0)) { // defect mode in case of clean site prob= ran4(&idum); stick_flag = 0; // no adsorption so far if (defects[x][y]==1) { // is the current site a defect site? if (prob > pdes_def) { // direct adsorption on clean defect site nstick++; x0=x; y0=y; stick_flag = 1; } else ndes++; // desorption from clean defect site }; if (defects[x][y]==0) { // is the current site a pristine site? if (prob > pdes_clean) { // direct adsorption on clean pristine site nstick++; x0=x; y0=y; stick_flag = 1; } else ndes++; // desorption from clean pristine site }; }; //********** extrinsic precursor in the case of a defect site if ( (def_label==0) & (lattice[x][y]==1) ) { // defect mode in case of a filled site prob= ran4(&idum); extr_flag=0; // no trapping in extrinsic precursor so far if ( (defects[x][y]==1) & (extr_flag==0) ) { // is the filled site a defect site? // i.e. defect mode in case of filled defect site extr_flag = 1; // yes, particle is now trapped in extrinsic precursor if (prob <= pdes_def_prec) { ndes++; } // desorption from filled defect site else { // trapping in extrinsic precursor l=0; hop_numbers= max_hops; while ( (lattice[x][y]==1) & (l pdes_clean) { // adsorption on clean site nstick++; x0=x; y0=y; stick_flag = 1; } else ndes++; // desorption from clean site }; }; //********** extrinsic precursor state in case of a site blocked by a preadsorbate if ((def_label==2)&((lattice[x][y]==1)|(defects[x][y]==1)) ) { //hit on site blocked by an adsorbate/preadsorbate prob= ran4(&idum); extr_flag=0; // adsorbate not trapped so far if ( (defects[x][y]==1) & (extr_flag==0) ) { // hit on preadsorbate site extr_flag = 1; if (prob <= pdes_def_prec) { ndes++; } // desorption else { //*** l=0; hop_numbers= max_hops; while (((lattice[x][y]==1)|(defects[x][y]==1))&(l pdes_clean) { // adsorption on clean site nstick++; x0=x; y0=y; stick_flag = 1; } else { if (intrin_hop>1) { hop_direction = (ran4(&idum)*8)+1; // intrinsic precursor hop make_a_hop(hop_direction, &x, &y); }; }; intrin_hop_number++; } /* while end */ if ( (intrin_hop_number > intrin_hop) & (stick_flag==0) ) ndes++; // extrinsic precursor loop start if ( lattice[x][y]==1) { // during the impact on filled site,the particle will be trapped in a // precursor state extr_flag = 1; l=0; hop_numbers= ran4(&idum)*max_hops; while ((lattice[x][y]==1) & (l pdes_occupied) { nstick++; x0=x; y0=y; } else { ndes++; // desorption from equilibrated state }; }; }; if ((extr_flag==0) & (stick_flag==0)) ndes++; *xx=x; *yy=y; *flag=stick_flag; } //********** main procedure start void monte_carlo() { float slope; /* initial slope of S(teta) vs. teta */ char answer[10]; int r_esponse; int i,j; /* rows, columns index */ int exposure; /* number of ettempts to adsorb a particle */ float prob; /* simulated desorption probability */ int x,y; /* lattice site */ int x0, y0; float temp; /* temporary */ int average_index; /* average results */ int temp_size; int hop_direction; int done; int ijk; float hop_numbers; float S_sum; int S_index; float cov_start, cov_end; /* for averaging */ float energy_start, energy_end; /* for averaging */ float A,B; float discuss; int discuss_flag; char filename[30]; int stick_flag; /* intrinsic precursor */ int extr_flag; /* particle was prappe in extrin. precursor state */ int jj; /* counts intrinsic precursor hops */ int x_line[3]; /* graphic output */ int y_line[3]; int x_index, y_index; int listen_index; int nlm; int x_liste[100]; /* adsoption sites of the particles during one averaging loop */ int y_liste[100]; char buffer[100]; char buffer2[100]; float initial_S; int color; char response[80]; /* handling of numeric overflows */ int status2; int iii, jjj; /* for thermalization prozedure */ //**** expose preadsorbates/defects for (i=1; i<=nx+3; i++) for (j=1;j<=ny+3; j++) defects[i][j]=0; if (( (def_label==0)| (def_label==2) ) & (def_cov > 0) ) pre_expose(); if (thermalize_defects_label==0) thermalize_defects(); //********** initialize variables for (i=1; i<=nx+4; i++) for(j=1; j<=ny+4; j++) lattice[i][j]=0; for (i=0;i<999;i++) { Stick[i]=99; coverage[i]=0; } lifetime = 0; curr_energy=0; // energy of the 1st (old) configuration new_energy=1; // energy of the next configuration index=1; // storing the averaged results average_index=0; // for averaging sum=0; // number of adsorbed particles S_sum=0; S_index=0; ndes=0; nstick=0; listen_index=0; // graphics output of the configuration //********** start monte carlo for (i=1; ( (sum < (nx*ny*saturation_coverage)) & (energy_loopndes) ); i++) { // EXPOSURE LOOP energy_loop=0; ndes=0; while ( ( new_energy >= curr_energy) & (energy_loop ]0,1[ */ y = (int)((ran4(&idum)*ny)+3); /* make sure that [3,ny+2] */ if (def_label==1) precursor(&x, &y, &stick_flag); if (def_label==0) defects_mode(&x, &y, &stick_flag); if (def_label==2) preadsorbates_mode(&x, &y, &stick_flag); }; //**********kinetic part start, energy constrained lattice[x][y]=1; // try to adsorb the particle, this defines the test configuration // backup of current configuration in the case of thermalization if (thermalize_label==1) { for (iii=1; iii<=(nx+2); iii++) { for (jjj=1; jjj<=(ny+2); jjj++) { temp_lattice[iii][jjj]=lattice[iii][jjj]; }; }; }; if (thermalize_label==1) thermalize(); if (boundary_flag==1) periodic_boundary (); if (grid_label==1) calc_energy(); // square lattice else calc_energy_hex(); // hexagonal lattice if ( (new_energy >= curr_energy) & (stick_flag==1) ) { // test configuration was not successful // keep current configuration if (thermalize_label==1) { for (iii=1; iii<=(nx+2); iii++) { for (jjj=1; jjj<=(ny+2); jjj++) { lattice[iii][jjj]=temp_lattice[iii][jjj]; }; }; }; lattice[x][y]=0; // start new energy cycle if (boundary_flag==1) periodic_boundary (); }; if ( (new_energy < curr_energy) ) { // test configuration was successful S_index++; A=(float)(nstick); B=(float)(ndes); S_sum = A/(A+B)+S_sum; if (S_index == 1) { cov_start = (float)(sum)/(float)((nx*ny)); energy_start = new_energy; } if (S_index == cov_average) { cov_end = (float)(sum)/(float)((nx*ny)); energy_end = new_energy; coverage[index] = (cov_end-cov_start)/2.0+cov_start; Stick[index] = S_sum/(float)(cov_average); convig_energy[index] = (energy_start/(sum+1-cov_average)+energy_end/sum)/2.0*(-1.0); if (index == 1) initial_S = Stick[index]; if (stick_norm != -1) Stick[index] = Stick[index]/initial_S*stick_norm; // print E vs. coverage printf("Theta = %f S=%f\n",coverage[index],Stick[index]); // print S vs. coverage if (index == 1) { coverage[0] = coverage[1]; Stick[0] = Stick[1]; }; // print configuration // print thermalized configuration S_index=0; index++; S_sum=0; listen_index=0; }; // coverage averaging, if }; // successful energy if }; // energy loop temp=curr_energy; curr_energy=new_energy; new_energy=temp; }; // exposure loop averaged_lifetime = (float)(lifetime)/(float)(i); // display configuration } int main(void) {monte_carlo();} **************************************************************************************************** //////////////////////////////////////////////////////////////////////////////////////////////////// **************************************************************************************************** Listing 4: //**********KMCS Langmuirian adsorption dynamics ste/bgh #include // standard "C" header #include // for simple mathematical functions #include // for simple input and output operations #include "nr.h" // from numerical recipes #define NRANSI // some compiler require this statement, // which indicates that standard "C" code follows int idum=(-1); // for the random number generator, include always this line in the examples // use the right rnd function rnd1 or rnd4 - that does not make a difference // or use one of the home made rnd functions, see the appendix //********** includes the final result float Coverage[500]; float AdsorptionProbability[500]; float AdsProbAverage[500]; float Realtime[500]; //********** model parameter int nx=50; // size of the lattice int ny=50; double Flux=1.0; // ML/sec normalised to nx, ny double Pclean=0.1; // desorption probability from clean sites int Average=10; // number of averaging MCS loops float SaturationCoverage=0.8; //********** program parameter int Lattice[200][200]; // adsorption sites of the lattice ; = 1 site is occupied; 0= empty site double Deposition[200][200]; // deposition prob. double Desorption[200][200]; // desorption prob. double Sclean; // adsorption prob. from clean sites double Dep; // deposition probability double Omega; // total transition probability double dt; // time increment double K; // for the linear selection procedure double Time; // time according to the current iteration double Told; // time according to the previous iteration int CurrentAveLoop; // current number of the averaging loops int Sum; // number of adsorbed particles int DepositionFlag = 0; // =1 a particle has been adsorbed ; =0 no adsorption int index=0; // index for graphic output controll int SaturationCoverageInt; // number of particles according to the saturation coverage void CalculateOmegaLD(int i, int j) { Omega-= Desorption[i][j]; Dep -= Deposition[i][j]; if (Lattice[i][j] == 0) { Desorption[i][j] = 0.0; Deposition[i][j] = Sclean; } else { Desorption[i][j] = 1.0; Deposition[i][j] = 0.0; } Omega += Desorption[i][j]; Dep += Deposition[i][j]; } void InitLD() { int i,j ; Omega = 0.0; Dep = 0.0; for (i=1; i<=nx; i++) { for (j=1;j<=ny; j++) { Desorption[i][j] = 0.0; Deposition[i][j] = 0.0 ; CalculateOmegaLD(i,j); } } } void DepositionProcessLD() { int i,j; int Flag=1; if (K < (Dep * Flux / (nx * ny))) { for (i=1; ( (i<=nx) & (Flag!=0) ); i++) { for (j=1;( (j<=nx) & (Flag!=0)); j++) { K -= (Deposition[i][j]*Flux/(nx*ny)); if (K <= 0) { Flag = 0; if (Lattice[i][j] == 0) { Lattice[i][j] = 1; Sum++; CalculateOmegaLD(i, j); } } } // for j } // for i DepositionFlag = 1; }// if K < ... end } void MainLoopLD() { int i,j; // loop index int SumPreview = 0; // number of adsorbed particles in previous loop double CurrentCoverage; double CurrentAdsProb; Time = 0; Told = 0.0; Sum = 0; for (i = 1; i <= nx; i++) { for (j = 1; j <= ny; j++) { Lattice[i][j] = 0; } } for (i = 1; i <= 200; i++) AdsProbAverage[i] = 0; InitLD(); while (Sum <= SaturationCoverageInt) { Omega += Flux; dt = -log(1.0 - ran4(&idum)) / Omega; Time += dt; K = ran4(&idum) * Omega; DepositionProcessLD(); Omega -= Flux; if ( ((int)(100*Sum) % (nx*ny)) == 0) { if ((DepositionFlag == 1) & (Sum > SumPreview)) { DepositionFlag = 0; SumPreview = Sum; CurrentCoverage = (double)(Sum) / (nx * ny); CurrentAdsProb = (nx*ny/100)/(Time-Told); index++; Coverage[index]=CurrentCoverage; AdsorptionProbability[index]=CurrentAdsProb; AdsProbAverage[index]= CurrentAdsProb + AdsProbAverage[index]; Realtime[index]=Time; printf("Loop=%d Time=%f Coverage=%f S=%f \n", CurrentAveLoop,Realtime[index], Coverage[index],AdsProbAverage[index]); Told = Time; } // if DepositionFlag==1 } // if int... } // while } void KmcsLD() { int i; SaturationCoverageInt = (int)(SaturationCoverage*nx*ny); CurrentAveLoop = 1; while (CurrentAveLoop <= Average) { index = 0; Sclean = 1.0 - Pclean; MainLoopLD(); CurrentAveLoop++; } for (i=1; i // standard "C" header #include // for simple mathematical functions #include // for simple input and output operations #include "nr.h" // from numerical recipes #define NRANSI // some compiler require this statement, // which indicates that standard "C" code follows int idum=(-1); // for the random number generator, include always this line in the examples //********** includes the final results float Coverage[500]; float AdsorptionProbability[500]; float AdsProbAverage[500]; float Realtime[500]; //********** model parameter int nx=50; // size of the lattice int ny=50; double Flux=1.0; // ML/sec normalised to nx, ny double Pclean=0.1 ; // desorption probability from clean sites int Average=10; // number of averaging MCS loops float SaturationCoverage=0.8; //********** program parameter int Lattice[50001]; // list of the adsorption sites of the lattice ; = 1 site is occupied; 0= empty site double Deposition[50001]; // deposition prob. list double Desorption[50001]; // desorption prob. list double Sclean; // adsorption prob. from clean sites double Dep; // deposition probability double Omega; // total transition probability double dt; // time increment double K; // for the linear selection procedure double Time; // time according to the current iteration double Told; // time according to the previous iteration int CurrentAveLoop; // current number of the averaging loops int Sum; // number of adsorbed particles int DepositionFlag = 0; // =1 a particle has been adsorbed ; =0 no adsorption int index=0; // index for graphic output control int SaturationCoverageInt; // number of particles according to the saturation coverage void CalculateOmegaLD(int c) { Omega-= Desorption[c]; Dep -= Deposition[c]; if (Lattice[c] == 0) { Desorption[c] = 0.0; Deposition[c] = Sclean; } else { Desorption[c] = 1.0; Deposition[c] = 0.0; } Omega += Desorption[c]; Dep += Deposition[c]; } void InitLD() { int i; Time = 0; Told = 0.0; Sum = 0; for (i = 1; i <= nx*ny; i++) Lattice[i] = 0; for (i = 1; i <= 499; i++) AdsProbAverage[i] = 0; Omega = 0.0; Dep = 0.0; for (i = 1; i <= nx * ny; i++) { Desorption[i] = 0.0; Deposition[i] = 0.0; CalculateOmegaLD(i);} } void DepositionProcessLD() { long int I=1,X=1; if (K < (Dep * Flux / (nx * ny))) { while (I == 1) { K -= (Deposition[X]*Flux/(nx*ny)); if (K <= 0) { I = 0; if (Lattice[X] == 0) { Lattice[X] = 1; Sum++; CalculateOmegaLD(X);} } X++; } // while end DepositionFlag = 1; }// if K < ... end } void MainLoopLD() { int SumPreview = 0; // number of adsorbed particles in previous loop double CurrentCoverage; double CurrentAdsProb; InitLD(); while (Sum <= SaturationCoverageInt) { Omega += Flux; dt = -log(1.0 - ran4(&idum)) / Omega; Time += dt; K = ran4(&idum) * Omega; DepositionProcessLD(); Omega -= Flux; if ( ((int)(100*Sum) % (nx*ny)) == 0) { if ((DepositionFlag == 1) & (Sum > SumPreview)) { DepositionFlag = 0; SumPreview = Sum; CurrentCoverage = (double)(Sum) / (nx * ny); CurrentAdsProb = (nx*ny/100)/(Time-Told); index++; Coverage[index]=CurrentCoverage; AdsorptionProbability[index]=CurrentAdsProb; AdsProbAverage[index]= CurrentAdsProb + AdsProbAverage[index]; Realtime[index]=Time; printf("Loop=%d Time=%f Coverage=%f S=%f\n", CurrentAveLoop, Realtime[index],Coverage[index],AdsProbAverage[index]); Told = Time; } // if DepositionFlag==1 } // if int... } // while } void KMCS_LD() { int i; SaturationCoverageInt = (int)(SaturationCoverage*nx*ny); CurrentAveLoop = 1; while (CurrentAveLoop <= Average) { index = 0; Sclean = 1.0 - Pclean; MainLoopLD(); CurrentAveLoop++; } for (i=1; i // standard "C" header #include // for simple mathematical functions #include // for simple input and output operations #include "nr.h" // from numerical recipes #define NRANSI // some compiler require this statement, // which indicates that standard "C" code follows float Coverage[500]; float AdsorptionProbability[500]; float AdsProbAverage[500]; float Realtime[500]; //********** model parameter int nx=100; // size of the lattice int ny=100; double Flux=1.0; // ML/sec normalised to nx, ny double Poccu=0.7; // desorption prob. from occupied sites double Pclean=0.6; // desorption probability from clean sites double Pdiff=1.2; // diffusion probability int Average=5; // number of averaging MCS loops float SaturationCoverage=0.3; //********** program parameter int Lattice[50001]; // list of the adsorption sites of the lattice = 1 site is occupied; 0= empty site int Precursor[50001]; // precursor state list int Neighbors[50001]; // includes the number of nearest neighbors for every lattice site int XR1[50001]; // neighbor list right neighbors int XL1[50001]; // left int YL1[50001]; // top int YR1[50001]; // bottom double Diffusion[50001]; // list including the diffusion probability for each adsorption site double Deposition[50001]; // deposition prob. list double Desorption[50001]; // desorption prob. list int idum=(-1); // for the Random number generator double Soccu; // adsorption prob. from occupied sites double Sclean; // adsorption prob. from clean sites double Dep; // deposition probability double Omega; // total transition probability double dt; // time increment double K; // for the linear selection procedure double Time; // time according to the current iteration double Told; // time according to the previous iteration int CurrentAveLoop; // current number of the averaging loops int Sum; // number of adsorbed particles int DepositionFlag = 0; // =1 a particle has been adsorbed; =0 no adsorption int index=0; // index for graphic output control int SaturationCoverageInt; // number of particles according to the saturation coverage void CalculateOmega(int c) { Omega-= (Desorption[c] + Diffusion[c]); Dep -= Deposition[c]; Neighbors[c] = Precursor[XL1[c]] + Precursor[XR1[c]] + Precursor[YL1[c]] + Precursor[YR1[c]]; if (Precursor[c] == 0) { Desorption[c] = 0.0; Diffusion[c] = 0.0; if (Lattice[c] == 0) Deposition[c] = Sclean; else Deposition[c] = Soccu; } else { Desorption[c] = Poccu; Deposition[c] = 0.0; Diffusion[c] = (4.0 - Neighbors[c]) * Pdiff; } Omega += (Desorption[c] + Diffusion[c]); Dep += Deposition[c]; } void INIT() { int I; Omega = 0.0; Dep = 0.0; for (I = 1; I <= nx * ny; I++) { Desorption[I] = 0.0; Deposition[I] = 0.0; Diffusion[I] = 0.0; CalculateOmega(I); } } void PeriodicBoundary() { long int I; for (I = 1; I <= nx*ny; I++) { XR1[I] = I + 1; XL1[I] = I - 1; YL1[I] = I - nx; YR1[I] = I + nx; } for (I = 1; I <= nx; I++) { YL1[I] = nx * (ny - 1) + I; YR1[nx * (ny - 1) + I] = I; } for (I = 1; I <= ny; I++) { XL1[1 + nx * (I - 1)] = nx * I; XR1[nx * I] = nx * (I - 1) + 1; } } void UPDATE(int a, int b) { int I2,I1; I1 = XL1[a]; CalculateOmega(I1); I1 = XR1[a]; CalculateOmega(I1); I1 = YL1[a]; CalculateOmega(I1); I1 = YR1[a]; CalculateOmega(I1); if (b > 0) { if (b==1) I1 = XL1[a]; if (b==2) I1 = XR1[a]; if (b==3) I1 = YL1[a]; if (b==4) I1 = YR1[a]; I2 = XL1[I1]; CalculateOmega(I2); I2 = XR1[I1]; CalculateOmega(I2); I2 = YL1[I1]; CalculateOmega(I2); I2 = YR1[I1]; CalculateOmega(I2); } } void DepositionProcess() { long int I=1,X=1; if (K < (Dep * Flux / (nx * ny))) { while (I == 1) { K -= (Deposition[X]*Flux/(nx*ny)); if (K <= 0) { I = 0; if (Lattice[X] == 0) { Lattice[X] = 1; Sum++; CalculateOmega(X);} else { Precursor[X] = 1; CalculateOmega(X);} UPDATE(X,0); } X++; } // while end DepositionFlag = 1; }// if K < ... end } void ChooseProcess() { int I,X; K -= Flux; I = 1; X = 1; while (I == 1) { K -= Desorption[X]; if (K < 0) { I = 0; Precursor[X] = 0; UPDATE(X,0); } else { K -= Diffusion[X]; if (K < 0) { K = ran4(&idum) * (4 - Neighbors[X]); Precursor[X] = 0; if (Precursor[XL1[X]] == 0) { K -= 1; if (K <= 0) { if (Lattice[XL1[X]] == 0) { DepositionFlag = 1; Sum++; Lattice[XL1[X]] = 1; } else Precursor[XL1[X]] = 1; CalculateOmega(X); UPDATE(X,1); I = 0; } } if ((Precursor[XR1[X]] == 0) & (I == 1)) { K -= 1; if (K <= 0) { if (Lattice[XR1[X]] == 0) { DepositionFlag = 1; Sum++; Lattice[XR1[X]] = 1; } else Precursor[XR1[X]] = 1; CalculateOmega(X); UPDATE(X,2); I = 0; } } if ((Precursor[YL1[X]] == 0) & (I == 1)) { K -= 1; if (K <= 0) { if (Lattice[YL1[X]] == 0) { DepositionFlag = 1; Sum++; Lattice[YL1[X]] = 1; } else Precursor[YL1[X]] = 1; CalculateOmega(X); UPDATE(X,3); I = 0; } } if ((Precursor[YR1[X]] == 0) & (I == 1)) { K -= 1; if (K <= 0) { if (Lattice[YR1[X]] == 0) { DepositionFlag = 1; Sum++; Lattice[YR1[X]] = 1; } else Precursor[YR1[X]] = 1; CalculateOmega(X); UPDATE(X,4); I = 0; } else { printf(" Error \n"); } } } } X++; } } void MainLoop() { int i; // loop index int SumPreview = 0; // number of adsorbed particles in previous loop double CurrentCoverage; double CurrentAdsProb; Time = 0; Told = 0.0; Sum = 0; PeriodicBoundary(); for (i = 1; i <= nx*ny; i++) Precursor[i] = 0; for (i = 1; i <= nx*ny; i++) Lattice[i] = 0; for (i = 1; i <= 499; i++) AdsProbAverage[index] = 0; INIT(); while (Sum <= SaturationCoverageInt) { Omega += Flux; dt = -log(1.0 - ran4(&idum)) / Omega; Time += dt; K = ran4(&idum) * Omega; if (K < Flux) DepositionProcess(); else ChooseProcess(); Omega -= Flux; if ( ((int)(100*Sum) % (nx*ny)) == 0) { if ((DepositionFlag == 1) & (Sum > SumPreview)) { DepositionFlag = 0; SumPreview = Sum; CurrentCoverage = (double)(Sum) / (nx * ny); CurrentAdsProb = (nx*ny/100)/(Time-Told); index++; Coverage[index]=CurrentCoverage; AdsorptionProbability[index]=CurrentAdsProb; AdsProbAverage[index]= CurrentAdsProb + AdsProbAverage[index]; Realtime[index]=Time; Told = Time; } // if DepositionFlag==1 } // if int... } // while } void KMCS() { int i; char answer[1]; SaturationCoverageInt = (int)(SaturationCoverage*nx*ny); CurrentAveLoop = 1; while (CurrentAveLoop <= Average) { printf("Loop No. %d \n",CurrentAveLoop); index = 0; Sclean = 1.0 - Pclean; Soccu = 1.0 - Poccu; MainLoop(); for (i=1; ib?a:b); } void Deposit() { // random deposition model char answer[1]; int i; // loop index int I; // adsorption site for (i = 1; i <= XSIZE; i++) X[i] = 0; // initialize 1D lattice N = 0; // initialize counter for w[N] in ANALYSIS() for (i=1; i <= t; i++) { // deposit t particles I = (int)((ran4(&idum)*XSIZE)+1); // select site for deposition at random if (ModeFlag==1) X[I]++; // deposit particle else X[I] = maximum(X[XL1[I]],maximum(X[I]+1,X[XR1[I]])); // ballistic deposition rule if (i % ((int)(ML)*XSIZE) == 0) { // analyze the interface every "ML" N++; // update counter ANALYSIS(); // do the analysis printf("deposited particles= %d interface width =%f \n",i,w[N]); // this line might be removed } } } void main() { int i,j; // loop index char answer[1]; // waiting for a key stroke PBC(); // apply periodic boundary conditions for (i = 1; i <= Averaging; i++) { printf("loop No. = %d \n",i); Deposit(); // start deposition model; for (j = 1; j <= N; j++) {wAverage[j] = wAverage[j]+w[j];}; } printf("\nAveraging result\n"); scanf("%s",answer); for (j = 1; j <= N; j++) { wAverage[j] = wAverage[j]/Averaging; printf("deposited particles =%d interface width=%f \n",j*(int)(ML)*XSIZE,wAverage[j]); }; scanf("%s",answer); } **************************************************************************************************** //////////////////////////////////////////////////////////////////////////////////////////////////// **************************************************************************************************** Listing 8-C: // Standard C #include #include #include #include "nr.h" // from numerical recipes #define NRANSI //********** Model parameter int nx = 50; // grid size int ny = 50; // grid size double Flux=1e-3; // in ML/s double Temperature=90; // in K double Enn=0.2; // latteral interaction energy double Eb=0.3; // binding energy, heat of adsorption int MaxPercentCoverage = 31; //********** program, controll variables int idum=(-1); // for the random number generator #define K_ESC '\x1B' // ??? double Ni = 1.0e+13; // effective attempt frequency double kB = 0.0000861734; // Boltzmann constant double Rate[5]; // constant table of rate int Lattice[150][150]; // lattice table int Isite[8000]; // sites adatom int Jsite[8000]; int Neighbours[8000]; // number of neighbours for adatom double HopRate[8000]; // rate hop int AtomNum; // number of deposited particles int PercentCoverage; // current coverage long int RunCounter; // MCS loop counter double Omega; // total transition rate double DepositionRate; double Coverage; double Time; int Di[4] = { -1, 0, 1, 0 }; // Ähh? int Dj[4] = { 0, 1, 0, -1 }; void StatisticResults( int PercCover ) // save results { FILE *f; int n; char FileName[40]; // Name of the file which includes the configuration at const coverage double Etot; // total energy of the configuration sprintf(FileName,"%d",PercCover); strcat( FileName, "Config"); // adds numerical value of PercCover to sting "Config" f = fopen(FileName,"wt"); // open configuration data file fprintf(f, "%d\n", AtomNum ); // save number of deposited particles for ( n = 1; n <= AtomNum; n++ ) fprintf(f, "%4d %4d\n", Isite[n], Jsite[n] ); // positions of occupied sites fclose(f); Etot = 0; for ( n = 1; n <= AtomNum; n++ ) Etot += Neighbours[n] * Enn; // calculates otal energy of configuration } int Xcoor( int i ) // ????? { return ( ( i + nx ) % nx ); } int Ycoor( int j ) // ????? { return ( ( j + ny ) % ny ); } int IsAtom( int i, int j ) { return ( Lattice[Xcoor(i)][Ycoor(j)] ); } void UpdateConfiguration( int n ) { int d; Neighbours[n] = 0; for ( d = 0; d < 4; d++ ) if ( IsAtom(Isite[n]+Di[d],Jsite[n]+Dj[d]) > 0 ) Neighbours[n]++; Omega -= HopRate[n]; HopRate[n] = ( 4 - Neighbours[n] ) * Rate[Neighbours[n]]; Omega += HopRate[n]; } void DepositionProcess(void) { int i, j, k, d; do { i = (int)( ran4(&idum) * (double)nx ); j = (int)( ran4(&idum) * (double)ny ); } while ( IsAtom(i,j) > 0 ); AtomNum++; Isite[AtomNum] = i; Jsite[AtomNum] = j; HopRate[AtomNum] = 0.0; Lattice[i][j] = AtomNum; UpdateConfiguration(AtomNum); for ( d = 0; d < 4; d++ ) if ( ( k = IsAtom(i+Di[d],j+Dj[d]) ) > 0 ) UpdateConfiguration(k); } void HopProcess( int n, int direct ) { int ia, ja, k, d; ia = Isite[n]; ja = Jsite[n]; Lattice[ia][ja] = 0; Isite[n] = Xcoor( ia+Di[direct] ); Jsite[n] = Ycoor( ja+Dj[direct] ); Lattice[Isite[n]][Jsite[n]] = n; for ( d = 0; d < 4; d++ ) { if ( ( k = IsAtom(ia+Di[d],ja+Dj[d]) ) > 0 ) UpdateConfiguration(k); if ( ( k = IsAtom(Isite[n]+Di[d],Jsite[n]+Dj[d]) ) > 0 ) UpdateConfiguration(k); } } void ChooseProcess(void) { int n, d; double RV, DeltaTime; if ( ( RunCounter % 1000L ) == 0L ) { Omega = 0.0; for ( n = 1; n <= AtomNum; n++ ) { HopRate[n] = 0.0; UpdateConfiguration(n); } } DepositionRate = (double)(nx*ny-AtomNum) * Flux; DeltaTime = -log( 1 - ran4(&idum) ) / ( Omega + DepositionRate ); Time = Time + DeltaTime; RV = ran4(&idum) * ( Omega + DepositionRate ); if ( RV >= Omega ) { DepositionProcess(); Coverage = (double)(AtomNum) / (double)(nx*ny); } else { n = 0; while ( RV >= 0 && n < AtomNum ) RV -= HopRate[++n]; d = 4; while ( RV < 0 && d > 0 ) { d--; if ( IsAtom( Isite[n]+Di[d], Jsite[n]+Dj[d] ) == 0 ) RV += Rate[Neighbours[n]]; } if ( IsAtom( Isite[n]+Di[d], Jsite[n]+Dj[d] ) == 0 ) HopProcess( n, d ); } } void SetAtom( void ) { FILE *f; int n, i, j, d, k, Num, ok; if ( ( f = fopen("INI.ATO","rt") ) != NULL ) { fscanf( f, "%d\n", &Num ); for ( n = 1; n <= Num; n++ ) { fscanf( f, "%d %d\n", &i, &j ); ok = ( i >= 0 ) && ( i < nx ) && ( j >= 0 ) && ( j < ny ) && ( IsAtom(i,j) == 0 ); if ( ok ) { AtomNum++; Lattice[i][j] = AtomNum; Isite[AtomNum] = i; Jsite[AtomNum] = j; HopRate[AtomNum] = 0; UpdateConfiguration(AtomNum); for ( d = 0; d < 4; d++ ) if ( ( k = IsAtom(i+Di[d],j+Dj[d]) ) > 0 ) UpdateConfiguration(k); } } fclose(f); } } void Initialize(void) { int i, j; for ( i = 0; i < 4; i++ ) Rate[i] = exp( log(Ni) - ( Eb + i*Enn )/(kB*Temperature) ); Rate[4] = 0; for ( i = 0; i < nx; i++ ) for ( j = 0; j < ny; j++ ) Lattice[i][j] = 0; AtomNum = 0; Omega = 0.0; Time = 0.0; SetAtom(); Coverage = (double)(AtomNum) / (double)(nx*ny); PercentCoverage = (int)(Coverage * 100) + 1; RunCounter = 0L; } void KMCS() { double x = sqrt(2); printf("KMCS start \n"); Initialize(); do { ChooseProcess(); RunCounter++; if ( ( Coverage*100 ) >= (double)PercentCoverage ) { printf("Coverage = %d\n",PercentCoverage); StatisticResults( PercentCoverage ); PercentCoverage++; } } while (( PercentCoverage <= MaxPercentCoverage ) ); StatisticResults( PercentCoverage ); } void main() {KMCS();} **************************************************************************************************** //////////////////////////////////////////////////////////////////////////////////////////////////// **************************************************************************************************** Listing 8-C & DOS: // DOS version #include // standard "C" header #include // for simple mathematical functions #include // for simple input and output operations #include "nr.h" // from numerical recipes #define NRANSI // some compiler require this statement, // which indicates that standard "C" code follows int nx, ny; // size of the grid Lattice, max 150x150 Rogo double Eb; // activation energy in eV double Enn; // interaction energy in eV double Flux; // flux in ML/s double Temperature; // temperature in Kelvin //* ***************** PROGRAM VARIABLES int Lattice[150][150]; // lattice table, 0=empty, 0 <> occupied int Isite[8000]; // x-coordinate of adatoms int Jsite[8000]; // y-coordinate of adatoms int Neighbors[8000]; // number of neighbors of adatom double HopRate[8000]; // diffusion rates int AtomNum; // number of adatoms deposited, i.e. adsorbed long int RunCounter; // number of times main loop passed double Sigma; // total diffusion rate = sum of hop rates of adatoms double DepositionRate; // number of adatoms hitting surface per second double Coverage; // fraction of monolayer covered = AtomNum / (nx * ny) double Time; // time elapsed since beginning of the simulation int PercentCoverage; // coverage multiplied by 100 and rounded to integral value // ***************** PROGRAM PARAMS & PRESETS double Ni = 1.0e+13; // effective attempt frequency double kB = 0.0000861734; // Boltzmann constant double Rate[5]; // table of rates for 0..4 neighbors int MaxPercentCoverage; // PercentCoverage by which program stops simulation int Di[4] = { -1, 0, 1, 0 }; // relative x-coordinate of surrounding sites int Dj[4] = { 0, 1, 0, -1 }; // relative y-coordinate of surrounding sites; these sites are numbered void StatisticResults( int PercCover ) { FILE *f; int n; char FileName[40]; double Etot; itoa( PercCover, FileName, 10 ); strcat( FileName, "KMC.ATO"); f = fopen(FileName,"wt"); for ( n = 1; n <= AtomNum; n++ ) fprintf(f, "%4d %4d\n", Isite[n], Jsite[n] ); fclose(f); Etot = 0; for ( n = 1; n <= AtomNum; n++ ) Etot += Neighbors[n] * Enn; f = fopen("KMC.ENE","at"); fprintf(f, "%15.8lf %6d %15.6le\n", Coverage, AtomNum, Etot/2 ); fclose(f); } void Results(void) { FILE *f; f = fopen("KMC.RES","wt"); fprintf(f, "lattice size nx = %6d\n", nx ); fprintf(f, " ny = %6d\n", ny ); fprintf(f, "adatom number = %6d\n", AtomNum ); fprintf(f, "coverage = %15.8lf\n", Coverage ); fprintf(f, "flux = %15.6le\n", Flux ); fprintf(f, "barrier energy Eb = %10.4lf\n", Eb ); fprintf(f, "lateral energy Enn = %10.4lf\n", Enn ); fprintf(f, "temperature = %10.2lf\n", Temperature ); fprintf(f, "process counter = %10ld\n", RunCounter ); fprintf(f, "time = %15.6le\n", Time ); fclose(f); } void SimulationProgress(void) { //* display key informations to the terminal screen RunCounter++; gotoxy( 40, 10 ); cprintf( "%10ld", RunCounter ); gotoxy( 40, 11 ); cprintf( "%6d", AtomNum ); gotoxy( 40, 12 ); cprintf( "%15lf", Coverage ); gotoxy( 40, 13 ); cprintf( "%15.6le", Time ); } main() { double x = sqrt(2); DeleteFiles(); Initialize(); do { ChooseProcess(); SimulationProgress(); if ( ( Coverage*100 ) >= (double)PercentCoverage ) { StatisticResults( PercentCoverage ); PercentCoverage++; } } while ( ( PercentCoverage <= MaxPercentCoverage ) && ( ! EscKeyPressed() ) ); StatisticResults( PercentCoverage ); Results();} void Initialize(void) { int i, j; ScanParam(); for ( i = 0; i < 4; i++ ) Rate[i] = exp( log(Ni) - ( Eb + i*Enn )/(kB*Temperature) ); Rate[4] = 0; for ( i = 0; i < nx; i++ ) for ( j = 0; j < ny; j++ ) Lattice[i][j] = 0; AtomNum = 0; Sigma = 0.0; Time = 0.0; SetAtom(); Coverage = (double)(AtomNum) / (double)(nx*ny); PercentCoverage = (int)(Coverage * 100) + 1; MaxPercentCoverage = 20; RunCounter = 0L; clrscr(); // DOS command which clears the screen gotoxy( 20, 10 ); cputs("Process counter = "); gotoxy( 20, 11 ); cputs(" Atom number = "); gotoxy( 20, 12 ); cputs(" Coverage = "); gotoxy( 20, 13 ); cputs(" Time = "); gotoxy( 20, 15 ); cputs("Break ... Esc"); } void ChooseProcess(void) { int n, d; double K, DeltaTime; if ( ( RunCounter % 1000L ) == 0L ) { // initialization Sigma = 0.0; for ( n = 1; n <= AtomNum; n++ ) { HopRate[n] = 0.0; UpdateConfiguration(n); } } DepositionRate = (double)(nx*ny) * Flux; DeltaTime = -log( 1 - Random() ) / ( Sigma + DepositionRate ); // do not forget the implement a random Time = Time + DeltaTime; // number generator assigned to the Random() function K = Random() * ( Sigma + DepositionRate ); if ( K >= Sigma ) { // handling deposition processes DepositionProcess(); Coverage = (double)(AtomNum) / (double)(nx*ny); } else { // handling diffusion events n = 0; while ( K >= 0 && n < AtomNum ) K -= HopRate[++n]; // first while loop d = 4; while ( K < 0 && d > 0 ) { // second while loop d--; if ( IsAtom( Isite[n]+Di[d], Jsite[n]+Dj[d] ) == 0 ) K += Rate[Neighbors[n]]; } if ( IsAtom( Isite[n]+Di[d], Jsite[n]+Dj[d] ) == 0 ) HopProcess( n, d ); } } void DepositionProcess(void) { int i, j, k, d; do { i = (int)( Random() * (double)nx ); j = (int)( Random() * (double)ny ); } while ( IsAtom(i,j) > 0 ); AtomNum++; Isite[AtomNum] = i; Jsite[AtomNum] = j; HopRate[AtomNum] = 0.0; Lattice[i][j] = AtomNum; UpdateConfiguration(AtomNum); for ( d = 0; d < 4; d++ ) if ( ( k = IsAtom(i+Di[d],j+Dj[d]) ) > 0 ) UpdateConfiguration(k); } int Xcoor( int i ) // return x-coordinate of adatom securing periodic boundary condition { return ( ( i + nx ) % nx ); } int Ycoor( int j ) // return y-coordinate of adatom securing periodic boundary condition { return ( ( j + ny ) % ny ); } void HopProcess( int n, int direct ) { int ia, ja, k, d; ia = Isite[n]; ja = Jsite[n]; Lattice[ia][ja] = 0; Isite[n] = Xcoor( ia+Di[direct] ); Jsite[n] = Ycoor( ja+Dj[direct] ); Lattice[Isite[n]][Jsite[n]] = n; for ( d = 0; d < 4; d++ ) { if ( ( k = IsAtom(ia+Di[d],ja+Dj[d]) ) > 0 ) UpdateConfiguration(k); if ( ( k = IsAtom(Isite[n]+Di[d],Jsite[n]+Dj[d]) ) > 0 ) UpdateConfiguration(k); } } void UpdateConfiguration( int n ) { int d; Neighbors[n] = 0; for ( d = 0; d < 4; d++ ) if ( IsAtom(Isite[n]+Di[d],Jsite[n]+Dj[d]) > 0 ) Neighbors[n]++; Sigma -= HopRate[n]; HopRate[n] = ( 4 - Neighbors[n] ) * Rate[Neighbors[n]]; Sigma += HopRate[n]; } **************************************************************************************************** //////////////////////////////////////////////////////////////////////////////////////////////////// **************************************************************************************************** Listing 9-FORTRAN & comments: program sigmaluca_bur_rev c corrected version washboard model c DECLARATIONS AND DIMENSIONING OF VARIABLES c declaration of variables double precision tetamin,tetamax double precision x10,y10,z10 double precision x1,y1,z1 double precision x20,y20,z20 double precision x2,y2,z2 double precision x3,y3,z3,x4,y4,z4,v3x,v3y,v3z,v4x,v4y,v4z double precision v1x0,v1y0,v1z0 double precision v1x,v1y,v1z double precision v2x0,v2y0,v2z0 double precision v2x,v2y,v2z double precision x1in,y1in,z1in,dt,en,en2,teta double precision x0,y0,z0 double precision x,y,z double precision diffen,enin,enfin,en10,en20,en1norm,en2norm double precision avexedes,avexedesn,aveoxdes,aveoxdesn double precision tetaout character*20 luca1,luca2 integer ndes,nrim,nx,ny,xeins,oxins,flad,ndir,nteta,index integer nocoll,nlost double precision sigmadir,deltasigmadir double precision vx0,vy0,vz0 double precision vx,vy,vz double precision xdes,ydes,deltasigmades double precision emin,emax,sigmades,ubin,x2in,y2in,z2in,zmin integer nr,flag,i1ag,j1ag,i2ag,j2ag,npart,nen double precision aveox,avexe,avesurf,avesurfdes double precision t,t12,t1ag,t2ag,tempo double precision en1x,en1y,en1z double precision r1,r2,r,a,sep double precision m1,m2,m,corr c dimensioning of matrices dimension aveox(20,20),avexe(20,20),avesurf(20,20) dimension avexedes(20,20),aveoxdes(20,20),avesurfdes(20,20) dimension avexedesn(20,20),aveoxdesn(20,20) dimension en(20,20),teta(20,20) dimension sigmades(20,20),sigmadir(20,20) dimension deltasigmades(20,20),deltasigmadir(20,20) dimension xdes(1000),ydes(1000) dimension x(16,16),y(16,16),z(16,16) dimension x0(16,16),y0(16,16),z0(16,16) dimension vx0(16,16),vy0(16,16),vz0(16,16) dimension vx(16,16),vy(16,16),vz(16,16) dimension tetaout(16,16,1000) c common statement for variables used in subroutines common/dati/ x10,y10,z10,x20,y20,z20,x0,y0,z0,v1x0,v1y0,v1z0, % v2x0,v2y0,v2z0,vx0,vy0,vz0,m1,m2,m,r1,r2,r c xe c setting of the arguments of random number generator function nx=16523 ny=98223 c initialize positions C INPUT OF THE PARAMETERS OF THE MODEL c opening file with input parameters open(unit=2,file='simuloluca_bur.dat',status='old') c min translational energy, max energy, number of energies read(2,*)emin,emax,nen c min angle of impingement, max angle, number of angles read(2,*)tetamin,tetamax,nteta c number of particles to be calculated for each enegry and angle read(2,*)npart c mass of the substrate atom, corugation parameter, separation c parameter c c m is an effective mass (in amu), which is larger than the mass of c substrate atoms. The larger the value of m, the lower the enegry c transfer to the substrate c c corr accounts for the surface corrugation and sep for the c distance between the adsorbed molecule and the surface c atoms, expressed in agstroms (10e-10 m) (see notes) read(2,*)m,corr,sep c binding energy in eV/molecule read(2,*)ubin c name of the output file containing the uotput fo the simulation read(2,5)luca1 c mass of the the projectile (m1=131 for Xe) and of the target c molecules (m2=32 for O2) (in amu) read(2,*)m1,m2 c radius in agstroms of projectile (r1=2.2 for Xe) and target (1.125 c for O2) read(2,*)r1,r2 c width of the square of the surface unit cell (2.88 A for Ag(001), c for which the lattice spacing is 2.88 * sqrt(2)=4.08 A) read(2,*)a c parameter specifying the surface mesh c 100 for (100) surfaces c 111 for (111) surfaces c 1001 for short bridge site (projectile velocity parallel to O-O c axis) c 1002 for short bridge site (projectile velocity normal to O-O c axis read(2,*)index c other output file with more details read(2,5)luca2 5 format(A) close(unit=2) c GEOMETRY OF THE SURFACE c calculate the radius of the surface atoms needed to simulate the c the corrugation c the surface is described as an array of balls of radius r, in c fixed positions which partly compenetrates each others c giving a washboard where the centers of the balls are a square c lattice of width a (for 100) surfaces, with the upmost point of c the sphere being higher by corr than the level at which balls c intersect each other (corr=0 is not allowd, but a c very small value (<< a) yiels a very large r, simulating thus c a flat surface ) r=0.5*((a/2)**2/corr+corr) c CONEVRTING INTO SI UNITS c initialize masses m1=m1*1.67e-27 m2=m2*1.67e-27 m=m*1.67e-27 r=r*1e-10 sep=sep*1e-10 c initialize dimensions r1=r1*1e-10 r2=r2*1e-10 c set the lateral dimension of the surface to 16X16 atoms nr=16 a=a*1e-10 c DEFINE POSITIONS OF THE ADSORBED MOLECULE x2in,y2in,z2in c they are put in the center of the lattice ... x2in=0 y2in=0 c at an height z2in from the z=0 in such a c way that the minimum separation between the hard walls c of the molecule and the surface atoms is sep c for (100) it is the cathet of the triangle formed c by r+r2+sep (hypothenusa) and a/2*(2**0.5) i.e. the half diagonal c of the surface unit cell c for fourfold hollow site on a square lattice if (index.eq.100) then z2in=((r+r2+sep)**2-(a/2*(2**0.5))**2)**0.5-r endif c for (111) surface the geometry is more complicated (TO BE CHECKED) c for threefold site on a 111 surface if (index.eq.111) then z2in=((r+r2+sep)**2-(a/4*(3**0.5))**2)**0.5-r endif c for the bridge sites the molecule is in between two atoms c and the calculation is simpler c for short bridge site Xe parallelo o-o caso b) c if (index.eq.1001) then z2in=((r+r2+sep)**2-(a/2)**2)**0.5-r endif c the same hold for the other direction of impingement c as the posiiton of the admolecule is the same c for short bridge site Xe ortogonale o-o caso c) if (index.eq.1002) then z2in=((r+r2+sep)**2-(a/2)**2)**0.5-r endif c OPEN OF THE FIRST OUTPUT FILE luca1 open(unit=3,file=luca1,status='new') c general informations about the simulation to check the values c used in the simulation and for memorandum when reading the c output file write(3,*)'program sigmaluca_bur_rev.for hard cube model' write(3,*)'During the simulation some particles may be' write(3,*)'lost as they go out of the lattice' write(3,*)'if such effect is larger than 2 % an ALERT appears' write(3,*)'m1,m2,m (amu)=',m1/1.67e-27,m2/1.67e-27,m/1.67e-27 write(3,*)'r1,r2,r,sep=(ang)',r1/1e-10,r2/1e-10,r/1e-10,sep/1e-10 write(3,*)'corr=',corr write(3,*)'zox=(ang)',z2in/1e-10 write(3,*)'nr=',nr write(3,*)'npart=',npart write(3,*) c from here on we have the header of the output file which ocntains c the order with which the outputs of the simulation are printed c c index of the simulation, energy, angle write(3,*)'JJ LL ENERGY TETA' c cross section for desorption and error write(3,*)' SIGMA DSIGMA' c sigma dir and error write(3,*)' SIGMA-DIR DSIGMA-DIR' c average values of Xe, oxyen, and surface energies write(3,*)'avexe(jj,ll),aveox(jj,ll),avesurf(jj,ll)' c average energy of Xe, oxygen and surface for desorption events write(3,*)'avexedes(jj,ll),aveoxdes(jj,ll),avesurfdes(jj,ll)' c average normal energy for Xe and oxygen for desorption events write(3,*)'avexedesn(jj,ll),aveoxdesn(jj,ll)' write(3,*) c OUTPUT TO THE SCREEN OF THE PARAMETERS USED IN THE c SIMULATION write(6,*)'m1,m2,m=',m1/1.67e-27,m2/1.67e-27,m/1.67e-27 write(6,*)'r1,r2,r,sep=',r1/1e-10,r2/1e-10,r/1e-10,sep/1e-10 write(6,*)'corr=',corr write(6,*)'zox=(ang)',z2in/1e-10 write(6,*)'nr=',nr write(6,*)'npart=',npart write(6,*) c OUTPUT ON SCREEN OF SIMULATION HEADERS AND ORDER OF OUTPUTS write(6,*)'JJ LL ENERGY TETA' write(6,*)' SIGMA DSIGMA' write(6,*)' SIGMA-DIR DSIGMA-DIR' write(6,*)'avexe(jj,ll),aveox(jj,ll),avesurf(jj,ll)' write(6,*)'avexedes(jj,ll),aveoxdes(jj,ll),avesurfdes(jj,ll)' write(6,*)'avexedesn(jj,ll),aveoxdesn(jj,ll)' write(6,*) c OPEN OF THE SECOND OUTPUT FILE luca2 open(unit=10,file=luca2,status='new') c summary of parameters employed for memorandum (meaning of symbols c is the same as before) write(10,*)'program sigmaluca_bur_rev.for hard cube model' write(10,*)'During the simulation some particles may be' write(10,*)'lost as they go out of the lattice' write(10,*)'if such effect is larger than 2 % an ALERT appears' write(10,*)'m1,m2,m=',m1/1.67e-27,m2/1.67e-27,m/1.67e-27 write(10,*)'r1,r2,r,sep=',r1/1e-10,r2/1e-10,r/1e-10,sep/1e-10 write(10,*)'corr=',corr c initial position of the target oxygen molecules write(10,*)'zox=(ang)',z2in/1e-10 c size of the lattice write(10,*)'nr=',nr c number of particles write(10,*)'npart=',npart write(10,*) c BEGINNING OF THE SIMULATION c cycle on nen different values of energy do 20000 jj=1, nen c cycle on nteta different values of the nagle do 15000 ll=1, nteta c if one energy only is wanted set the enegry to the minimum value if (nen.eq.1) then en(1,ll)=emin else c otherwise calculate actual energy en(jj,ll)=emin+(emax-emin)*(jj-1)/(nen-1) endif c if one angel is wanted set it to the minimum value if (nteta.eq.1) then teta(jj,1)=tetamin else c otherwise calculate actual angle value teta(jj,ll)=tetamin+(tetamax-tetamin)*(ll-1)/(nteta-1) endif c initialize coutning parameters c numebr of desorbing particles ndes=0 c CHECK ndir=0 c number of remaining particles nrim=0 c number of particle which do not collide nocoll=0 c numebr of particles lost nlost=0 c number of Xe (xeins) and oxygen (oxins) which for some errore go c inside the surface (z<0) xeins=0 oxins=0 c initialize matrices containing the uotput of the simulation c cross section for desorption sigmades(jj,ll)=0 c sigmadir CHECK sigmadir(jj,ll)=0 c average ox energy aveox(jj,ll)=0 c average surface energy for desorpiton events, the same for ox, Xe avesurfdes(jj,ll)=0 aveoxdes(jj,ll)=0 avexedes(jj,ll)=0 c as above but for normal energy aveoxdesn(jj,ll)=0 avexedesn(jj,ll)=0 avexe(jj,ll)=0 avesurf(jj,ll)=0 c cycle on number of particles do 10000 ii=1, npart c error check on total number of desorbing and remaining particles c if it is wrong error message and go out of the cycle if ((ndes+nrim).gt.(npart)) then write(6,*)'ERROR on total particles numbers' write(3,*)'ERROR on total particles numbers' write(10,*)'ERROR on total particles numbers' goto 10010 endif c c initialize energy of the incoming Xe atom along x,y,z en1y=0 en1x=en(jj,ll)*(sin(3.1415*teta(jj,ll)/180.))**2 en1z=en(jj,ll)*(cos(3.1415*teta(jj,ll)/180.))**2 c set the initial z coordinate of Xe (z10) z10=100 c generate x (x10) and y (y10) impact coordinates for Xe c over a square of width 20 angstroms centered in (0,0) c NOTE THAT IF THE SQUARE IS CHANGED ALSO THE NORMALIZATION c factor 400 used in calculating sigma at line 10010 MUST c be MODIFIED x10=-10+20*ran(nx) y10=-10+20*ran(ny) c intialize x1in,y1in x1in=x10 y1in=y10 c translate the initial position in such a way that the impact c coordinate is at the desired position (necessary for non normal c incidence) x10=x10-z10*tan(3.1415*teta(jj,ll)/180) c convert into angstroms x10=x10*1e-10 y10=y10*1e-10 z10=z10*1e-10 z1in=z10 c initialize x1,y1,z1 (Xe next time postition) x1=x10 y1=y10 z1=z10 c set oxygen coordinates (x20,y20,z20) to the initial adsorbed c positions (x2in,y2in,z2in) x20=x2in y20=y2in z20=z2in c initialize x2,y2,z2 (oxygen next time position) x2=x20 y2=y20 z2=z20 c BUILD THE SQUARE LATTICE (z=0 corresponds to the upper point c of the surface ball, c i.e. the plane tangent to the surface has coordinate z=0) c silver atoms c cicle on x and y coordinates (indexes i,j respectively) do 10 i=1, nr do 9 j=1, nr if (index.eq.100) then c fourfold hollow x0(i,j)=-0.5*(nr-1)*a+(i-1)*a y0(i,j)=-0.5*(nr-1)*a+(j-1)*a z0(i,j)=-r endif c top position c x0(i,j)=0 c y0(i,j)=0 c z0(i,j)=-r c the same for (111) lattice if (index.eq.111) then c threefold site x0(i,j)=-0.5*a+(i-1)*a+(j-1)*a/2-(nr*0.5-1)*a-(nr*0.5-1)*a/2 y0(i,j)=a*(3.**0.5)/4+(j-1)*a*(3.**0.5)/2- % (nr*0.5-1)*a*(3.**0.5)/2 z0(i,j)=-r endif c ... and for bridge sites c short bridge site xe parallelo o-o caso b) if (index.eq.1001) then c short bridge x0(i,j)=-0.5*(nr-1)*a+(i-1)*a y0(i,j)=-0.5*(nr-1)*a+(j-1)*a+0.5*a z0(i,j)=-r endif c short bridge site xe ortogonale o-o caso c) if (index.eq.1002) then c short bridge x0(i,j)=-0.5*(nr-1)*a+(i-1)*a+0.5*a y0(i,j)=-0.5*(nr-1)*a+(j-1)*a z0(i,j)=-r endif c initialzie matrices also (x,y,z) silver next time positions x(i,j)=x0(i,j) y(i,j)=y0(i,j) z(i,j)=z0(i,j) 9 continue 10 continue c initialize velocities c xe c converts energy to SI en1x=en1x*1.6e-19 en1y=en1y*1.6e-19 en1z=en1z*1.6e-19 c calculate intial velocities for Xe (v1x0,v1y0,v1z0) v1x0=(2*en1x/m1)**0.5 v1y0=(2*en1y/m1)**0.5 c v1z0 is negative as the atoms is directed toward the surface v1z0=-(2*en1z/m1)**0.5 c initialize next time velocities for Xe (v1x,v1y,v1z) v1x=v1x0 v1y=v1y0 v1z=v1z0 c and for Oxygen c initializie initial time velocities for Oxygen v2x0=0 v2y0=0 v2z0=0 c and next time velocities for oxygen v2x=v2x0 v2y=v2y0 v2z=v2z0 c initilaize initial and next time velocities for silver atoms c cicles on the x,y array do 20 i=1, nr do 19 j=1, nr vx0(i,j)=0 vy0(i,j)=0 vz0(i,j)=0 vx(i,j)=0 vy(i,j)=0 vz(i,j)=0 19 continue 20 continue c c INTEGRATE EQUATION OF MOTION c set initial time tempo to 0 tempo=0 c calculate initial oxygen (en20) and Xe (en10) energies en20=0.5*m2*(v2x0**2+v2y0**2+v2z0**2) en10=0.5*m1*(v1x0**2+v1y0**2+v1z0**2) c cccccccccccccc c integrate equation of motions tempo=0 c set control variable flad to 0 flad=0 c cicle on collisio events do 1000 kk=1, 300 c determine the first collision event c between xe and ag c intialize t of collision 1 (xe) - Ag 1ag to a very high value c corresponding to a time fo collision in the far future c this value must be much larger than c (initial Xe distance from the surface)/(mimimum initial velocity c employed in the simulations) c if the program is used under very different conditions c t1ag has to be changed to satisfy such condition t1ag=1000000 c set control variable flag to 0 flag=0 c initilaize dt=0 dt=0 c c calculate using subroutine tim the time of collision of Xe c with each silver atom do 100 i=1, nr do 99 j=1, nr call tim(x10,y10,z10,x0(i,j),y0(i,j),z0(i,j),v1x0,v1y0,v1z0, % vx0(i,j),vy0(i,j),vz0(i,j),r1,r,t) c determine at which time (t1ag) and with each atom (i1ag,j1ag) c the first xe-ag collision occurs if (t.lt.t1ag) then t1ag=t i1ag=i j1ag=j endif 99 continue 100 continue c do the same for collison of xe with oxygen c between xe and 2 c set t12 to a large value corresponding to a Xe-Ox collision in c the far future (see comments above for t1ag and modify c t12 accordignly if necessary) t12=1000000 c calculate collision between Xe and Ox call tim(x10,y10,z10,x20,y20,z20,v1x0,v1y0,v1z0,v2x0,v2y0,v2z0, % r1,r2,t) c set t12 to the calculated value ( a time t longer than t12 c woudl imply no collision) if (t.lt.t12) then t12=t endif c as above for collisions between ox (2) and Ag (t2ag) c between ox and ag c set t2ag to a high values (see comments above for t1ag and t12) t2ag=1000000 c search for the time t of collision of ox with ag on the lattice do 200 i=1, nr do 199 j=1, nr call tim(x20,y20,z20,x0(i,j),y0(i,j),z0(i,j),v2x0,v2y0,v2z0, % vx0(i,j),vy0(i,j),vz0(i,j),r2,r,t) c determine the time t2ag of the first ox-ag collision c and the atom (i2ag,j2ag) at which the first collision occurs if (t.lt.t2ag) then t2ag=t i2ag=i j2ag=j endif 199 continue 200 continue c DETERMINE WHICH TYPE OF COLLISION OCCURS FIRST c write(6,*)'t1ag,t12,t2ag',t1ag,t12,t2ag if ((t1ag.lt.t2ag).and.(t1ag.lt.t12)) then c write(6,*)'urto tra xe e ag',t1ag c set control variable flag to 1 (= Xe-Ag collision) flag=1 dt=t1ag c increase time tempo by dt=t1ag tempo=tempo+dt endif if ((t12.lt.t1ag).and.(t12.lt.t2ag)) then c write(6,*)'urto tra xe e ox',t12 c set contorl variavle flag to 2 (= Xe-Ox collision) flag=2 dt=t12 c increase time tempo by dt=t12 tempo=tempo+dt endif if ((t2ag.lt.t1ag).and.(t2ag.lt.t12)) then c write(6,*)'urto tra ox e ag',t2ag c set control varibale flag to 3 (= Ox -Ag collision) flag=3 dt=t2ag c increase tyme tempo by dt=t2ag tempo=tempo+dt endif if (kk.eq.1) then c for the first collision set control variable flad c equal to flag (it ocntains the nature of the first collision c and is 0 for no collision) flad=flag endif c eliminate collision occurring after an unreasonable time if ((dt.gt.10000).and.(kk.eq.1)) then flad=0 endif c write(6,*)'jj,ii,kk,flag,flad',jj,ii,kk,flag,flad c if no collsiion is occurring than set kk=299 which ensures that c at the end of the cyle the program go to other initial positions c increase the counting nocoll if (flag.eq.0) then dt=1e-10 kk=299 c write(6,*)'No more collision' endif c the cotnrol variable flad is set equal to flag at the first c collision, becoming 1, 2 or 3 depending on the collision type. c If it is still 0 at this tage it m eans that no collision has c occured this will ocntirbute to an ALERT later if the number of c such events is larger than 2 % of the total number if (flad.eq.0) then nocoll=nocoll+1 endif c IN CASE OF COLLISION EVENT INTEGRATION O EQUATION c calculate new positions x1=x10+v1x0*dt y1=y10+v1y0*dt z1=z10+v1z0*dt x2=x20+v2x0*dt y2=y20+v2y0*dt z2=z20+v2z0*dt c calculate new velocIties C IN CASE OF NO COLLISIONS mantain old velocities if (flag.eq.0) then v1x=v1x0 v1y=v1y0 v1z=v1z0 v2x=v2x0 v2y=v2y0 v2z=v2z0 c go and check if penetration inside the surface has occurred goto 290 endif c for type 1 collision CALCULATE NEW COORDINATES c AND VELOCITIES USIGN THE SUBROUTINE HARD(...) c The subsrutine reads initial position for xe (x10,y10,z10), for c the silver atom involved (x0(i1ag,j1ag),y0(),z0) and velocities c (v1x0,v1y0,v1z0) for Xe and the Ag involved and yields c new position (x3,y3,z3) and velocity (v3x,v3y,v3z) for Xe and for c ag atom involved (x4,y4,z4) and (vx4,vy4,vz4) if (flag.eq.1) then call hard(x10,y10,z10,x0(i1ag,j1ag),y0(i1ag,j1ag),z0(i1ag,j1ag), % v1x0,v1y0,v1z0,vx0(i1ag,j1ag),vy0(i1ag,j1ag),vz0(i1ag,j1ag), % x3,y3,z3,x4,y4,z4,v3x,v3y,v3z,v4x,v4y,v4z,m1,m,r1,r,flag) c set next time values of velocity equal to the computed values v1x=v3x v1y=v3y v1z=v3z v2x=v2x0 v2y=v2y0 v2z=v2z0 c check for penetration inside goto 290 endif c for type 2 collision CALCULATE NEW COORDINATES c AND VELOCITIES USIGN THE SUBROUTINE HARD(...) c meanign of symbols as for case 1, with *1 and *2 Xe and Ox c inputs, and *3 and *4 xe and ox outputs if (flag.eq.2) then call hard(x10,y10,z10,x20,y20,z20, % v1x0,v1y0,v1z0,v2x0,v2y0,v2z0, % x3,y3,z3,x4,y4,z4,v3x,v3y,v3z,v4x,v4y,v4z,m1,m2,r1,r2,flag) c set next time velocitis to the calculated values v1x=v3x v1y=v3y v1z=v3z v2x=v4x v2y=v4y v2z=v4z goto 290 endif c for type 3 collision CALCULATE NEW COORDINATES c AND VELOCITIES USIGN THE SUBROUTINE HARD(...) c The subsrutine reads initial position for ox (x10,y10,z10), for c the silver atom involved (x0(i1ag,j1ag),y0(),z0) and velocities c (v1x0,v1y0,v1z0) for ox and the Ag involved and yields c new position (x3,y3,z3) and velocity (v3x,v3y,v3z) for ox and for c ag atom involved (x4,y4,z4) and (vx4,vy4,vz4) if (flag.eq.3) then call hard(x20,y20,z20,x0(i2ag,j2ag),y0(i2ag,j2ag),z0(i2ag,j2ag), % v2x0,v2y0,v2z0,vx0(i2ag,j2ag),vy0(i2ag,j2ag),vz0(i2ag,j2ag), % x3,y3,z3,x4,y4,z4,v3x,v3y,v3z,v4x,v4y,v4z,m2,m,r2,r,flag) c set next time values of velocity equal to the computed values v1x=v1x0 v1y=v1y0 v1z=v1z0 v2x=v3x v2y=v3y v2z=v3z c check for penetration inside goto 290 endif c row 290 CONTROLS IF penetration HAS OCCURRED and going on with the c cycle 10000 on npart after increasing the number of xeins (if it c is Xe which is goes inside) or oxins (if penetration involves ox) c by 1 unit 290 if (z1.lt.0) then xeins=xeins+1 goto 10000 endif if (z2.lt.0) then oxins=oxins+1 goto 10000 endif c IF PENETRATION HAS NOT OCCURED IT searches for hte next collision c and reset initial (*0) condition to the present state value (*1) c to go on with the next collision event and so on c reset boundary conditions x10=x1 y10=y1 z10=z1 x20=x2 y20=y2 z20=z2 v1x0=v1x v1y0=v1y v1z0=v1z v2x0=v2x v2y0=v2y v2z0=v2z c write(6,*)'tempo',tempo c write(6,*)'coordinate of xe',x10/1e-10,y10/1e-10,z10/1e-10 c write(6,*)'coordinate of ox',x20/1e-10,y20/1e-10,z20/1e-10 v10=(v1x0**2+v1y0**2+v1z0**2)**0.5 v20=(v2x0**2+v2y0**2+v2z0**2)**0.5 c if velocity of Xe is zero then... if (v10.eq.0) then goto 900 else c write(6,*)'xe velocit dir',v1x0/v10,v1y0/v10,v1z0/v10 endif c check if ox too has stopped 900 if (v20.eq.0) then goto 902 else c write(6,*)'ox velocit dir',v2x0/v20,v2y0/v20,v2z0/v20 endif c if the z coordinate of Xe (z10) is larger than c the starting value (z1in) it means Xe is leaving the surface c and there is no need for further searches c Then the program goes out of the cycle 1000 searching for c collisions going to line 1001 902 if (z10.gt.z1in) then c write(6,*)'z10 > z1in ',z10,z1in goto 1001 endif c if in any case no colliision has occured do the same exiting the c cycle 1000 goin to line 1001 if (flad.eq.0) then nlost=nlost+1 write(3,*)'LOST PARTICLE' write(3,*)'1 coordinate ',x10,y10,z10 write(3,*)'1 velocity',v1x0,v1y0,v1z0 write(3,*)'2 coordinate ',x20,y20,z20 write(3,*)'2 velocity',v2x0,v2y0,v2z0 write(3,*)'t12,t1ag,t2ag,t',t12,t1ag,t2ag,t goto 1001 endif 1000 continue c end of cycle seacrhing for collisions c line 1001 reached after 300 collisions c or after proving that no collision has occured 1001 en2norm=0.5*m2*v2z**2 en1norm=0.5*m1*v1z**2 en1=0.5*m1*(v1x**2+v1y**2+v1z**2) en2=0.5*m2*(v2x**2+v2y**2+v2z**2) c increase summing variables (which were set to 0 c before starting the cycle 10000 for simulation) c to compute the average of c c oxygen final (suffix 2) energy aveox(jj,ll)=aveox(jj,ll)+en2 c xe final energy avexe(jj,ll)=avexe(jj,ll)+en1 c surface final energy avesurf(jj,ll)=avesurf(jj,ll)+(en10-en1-en2) c CONDITION FOR OXYGEN DESORPTION C NORMAL ENERGY OF OXYGEN (en2norm) HIGHER THAN THE BINDING ENERGY c (ubin) AND DIRECTION OF VELOCITY LEAVING THE SURFACE (v2z>0) if ((en2norm.gt.(ubin*1.6e-19)).and.(v2z.gt.0)) % then c increase the variable ndes which counts the desorbed particles by c 1 unit ndes=ndes+1 c check if that happend after xe -ox collision c in this case increase the variable ndir which counts direct c desorption events by 1 unit if (flad.eq.2) then ndir=ndir+1 endif c write(6,*)'desorbed' c write on the screen such event write(6,*)'ox velocit dir',v2x0/v20,v2y0/v20,v2z0/v20 c calculate outcoming angle tetaout and put in a matrix c where the first and second index are x,y, running indexes c on the square lattice and the third index run over ndes tetaout(jj,ll,ndes)=180./3.1415* % acos(v2z/(v2x**2+v2y**2+v2z**2)**0.5) c write it on the screen write(6,*)'tetaout=',v2z/(v2x**2+v2y**2+v2z**2)**0.5, % tetaout(jj,ll,ndes) c increase sum variable for desorbing events for ... c surface energy avesurfdes(jj,ll)=avesurfdes(jj,ll)+(en(jj,ll)*1.6e-19)-0.5*m2 % *(v2x**2+v2y**2+v2z**2)- % 0.5*m1*(v1x**2+v1y**2+v1z**2) c oxygen energy aveoxdes(jj,ll)=aveoxdes(jj,ll)+0.5*m2*(v2x**2+v2y**2+v2z**2) c normal oxygen energy aveoxdesn(jj,ll)=aveoxdesn(jj,ll)+0.5*m2*(v2z**2) c xe energy avexedes(jj,ll)=avexedes(jj,ll)+0.5*m1*(v1x**2+v1y**2+v1z**2) c normal xe enegry avexedesn(jj,ll)=avexedesn(jj,ll)+0.5*m1*(v1z**2) c take memory of initial conditions leading to c such events by putting the initial coordinates for Xe c x1in,y1in in a vector (xdes,ydes) whihc contains ndes elements xdes(ndes)=x1in ydes(ndes)=y1in else c otherwise count the oxygen as remaining on the surface and go on c by increasing the variable nrim whcih counts remaining molecules c by 1 unit nrim=nrim+1 c write(6,*)'not desorbed' endif 10000 continue c error check on total number of desorbing and remaining particles c if it is wrong error message and go out of the cycle if ((ndes+nrim).gt.(npart)) then write(6,*)'ERROR on total particles numbers (excess)' write(6,*)'ndes,nrim,ndes+nrim',ndes,nrim,ndes+nrim write(3,*)'ERROR on total particles numbers (excess)' write(3,*)'ndes,nrim,ndes+nrim',ndes,nrim,ndes+nrim write(10,*)'ERROR on total particles numbers (excess)' write(10,*)'ndes,nrim,ndes+nrim',ndes,nrim,ndes+nrim endif c error check on total number of desorbing, penetrated and c remaining particles c if it is wrong error message and go out of the cycle if ((ndes+nrim+oxins+nocoll+nlost).lt.(npart)) then write(6,*)'ERROR on total particles numbers (defect)' write(6,*)'ndes,nrim,oxins,TOT',ndes,nrim,oxins,ndes+nrim+oxins write(3,*)'ERROR on total particles numbers (defect)' write(3,*)'ndes,nrim,oxins,TOT',ndes,nrim,oxins,ndes+nrim+oxins write(10,*)'ERROR on total particles numbers (defect)' write(10,*)'ndes,nrim,oxins,TOT',ndes,nrim,oxins,ndes+nrim+oxins write(6,*)'particle no colliding=',nocoll write(3,*)'particle no colliding=',nocoll write(10,*)'particle no colliding=',nocoll write(6,*)'particle lost',nlost write(3,*)'particle lost',nlost write(10,*)'particle lost',nlos endif c end of the cycle on npart c CALCULATE THE CROSS SECTIONS c by normalizing the summing variables c the division by 400 arises form the fact that initial c conditions over a square of 20 X 20 angstroms were chosen c when initialising x1in and y1n at the beginning c the values are thus directly in angtroms**2 c if some penetration event has occurred or if some particle is lost c the program check if the total number of suspect events is c larger than 1 % and gives in this case an ALERT indicating c that some problems may exist with the model which need further c ATTENTION c Penetration in the bulk may indicate that there are holes between c the atoms which modle the surface c Lost particles may indicate that the set of initial parameters c is too wide for the dimension of the array of atoms so that c particle nether collide with it 10010 sigmades(jj,ll)=ndes*1.0/((ndes+nrim)*1.0)*400.0 sigmadir(jj,ll)=ndir*1.0/((ndes+nrim)*1.0)*400.0 if (((nocoll+oxins+xeins+nlost)*1.0).gt.(0.02*npart)) then write(6,*)'nocoll=',nocoll write(6,*)'oxins=',oxins write(6,*)'xeins=',xeins write(6,*)'ndes,nrim=',ndes,nrim write(6,*)'npart=',npart write(6,*)' ALERT ! TOO PARTICLES INSIDE OR LOST' write(6,*)' TAKE RESULT WITH CARE POSSIBLE FLAW !!' write(3,*)'nocoll=',nocoll write(3,*)'oxins=',oxins write(3,*)'xeins=',xeins write(3,*)'ndes,nrim=',ndes,nrim write(3,*)'nlost=',nlost write(3,*)'npart=',npart write(3,*)(nocoll+oxins+xeins+nlost)*1.0 write(3,*)(0.02*npart) write(3,*)' ALERT ! TOO PARTICLES INSIDE OR LOST' write(3,*)' TAKE RESULT WITH CARE POSSIBLE FLAW !!' write(10,*)'nocoll=',nocoll write(10,*)'oxins=',oxins write(10,*)'xeins=',xeins write(10,*)'ndes,nrim=',ndes,nrim write(10,*)'npart=',npart write(10,*)' ALERT ! TOO PARTICLES INSIDE OR LOST' write(10,*)' TAKE RESULT WITH CARE POSSIBLE FLAW !!' write(6,*)'Penetration in the bulk may indicate that there' write(6,*)'are holes between the atoms which model the surface' write(6,*)'Lost particles may indicate that the set of initial' write(6,*)'parameters is too wide for the dimension of the' write(6,*)'array of atoms so that some particle never reach it' write(3,*)'Penetration in the bulk may indicate that there' write(3,*)'are holes between the atoms which model the surface' write(3,*)'Lost particles may indicate that the set of initial' write(3,*)'parameters is too wide for the dimension of the' write(3,*)'array of atoms so that some particle never reach it' endif c calculate error, assumign a Poisson statistic deltasigmades(jj,ll)=(ndes**0.5)*1.0/((ndes+nrim)*1.0)*400.0 deltasigmadir(jj,ll)=(ndir**0.5)*1.0/((ndes+nrim)*1.0)*400.0 c devide sum variable by the number of particles c to obtain AVERAGE ENERGIES aveox(jj,ll)=aveox(jj,ll)/((ndes+nrim)*1.6e-19) avexe(jj,ll)=avexe(jj,ll)/((ndes+nrim)*1.6e-19) avesurf(jj,ll)=avesurf(jj,ll)/((ndes+nrim)*1.6e-19) c if some particle has desorbed compute also c average values for desorbing ensambles if (ndes.gt.0) then aveoxdes(jj,ll)=aveoxdes(jj,ll)/(ndes*1.0)/1.6e-19 aveoxdesn(jj,ll)=aveoxdesn(jj,ll)/(ndes*1.0)/1.6e-19 avexedes(jj,ll)=avexedes(jj,ll)/(ndes*1.0)/1.6e-19 avexedesn(jj,ll)=avexedesn(jj,ll)/(ndes*1.0)/1.6e-19 avesurfdes(jj,ll)=avesurfdes(jj,ll)/(ndes*1.0)/1.6e-19 else c otherwise SET THE CORRESPONDING VALUES TO 0 aveoxdes(jj,ll)=0 aveoxdesn(jj,ll)=0 avexedes(jj,ll)=0 avexedesn(jj,ll)=0 avesurfdes(jj,ll)=0 endif c output to file 10 write(10,*)jj,ll,en(jj,ll),teta(jj,ll),ndes c IF DESORBED PARTICLES EXIST, THEN CALCULATE DESORPTION c angle by the subsoutine conto acting on the vector (xdes,ydes) c of dimension ndes build above if (ndes.gt.0) then call conto(jj,ll,ndes,tetaout) else write(10,*) endif c WRITE OUTPUT TO SCREEN AND TO THE FIRST OUTPUT FILE write(10,*)'----------' write(6,*)jj,ll,en(jj,ll),teta(jj,ll) write(6,*)sigmades(jj,ll),deltasigmades(jj,ll) write(6,*)sigmadir(jj,ll),deltasigmadir(jj,ll) write(6,*)'---------------' write(3,*)jj,ll,en(jj,ll),teta(jj,ll) write(3,*)sigmades(jj,ll),deltasigmades(jj,ll) write(3,*)sigmadir(jj,ll),deltasigmadir(jj,ll) write(3,*)avexe(jj,ll),aveox(jj,ll),avesurf(jj,ll) write(3,*)avexedes(jj,ll),aveoxdes(jj,ll),avesurfdes(jj,ll) write(3,*)avexedesn(jj,ll),aveoxdesn(jj,ll) write(3,*) 15000 continue c end cicle on angle 20000 continue c end cicle on enegry close(unit=3) close(unit=10) c CLOSE OUTPUT FILES AND END OF THE MAIN PROGRAM stop end c ----------------- c SUBROUTINES c ----------------- c subroutine conto(i,j,n,angolo) integer i,j,n,nt double precision angolo,t,tmin,tmax dimension angolo(16,16,1000),nt(9) do 10 l=1, 9 nt(l)=0 10 continue do 100 k=1, n do 80 l=1, 9 tmin=(l-1)*10. tmax=l*10. c determine the cell within which the out nagle has to be put if ((angolo(i,j,k).gt.tmin).and.(angolo(i,j,k).lt.tmax)) then nt(l)=nt(l)+1 endif 80 continue 100 continue do 110 l=1, 9 write(10,*)(l-0.5)*10.,nt(l)/(n*1.0) 110 continue end c ------------- c the suboutine distrib is not called anywhere c but it can be useful for detailed analysis of the desorbing c particles subroutine distrib(i,j,nd,angolo) c writes in a file distribuzione.dat the matrix angolo c folding it into nd square matrices 10X10 double precision angolo integer i,j,nd dimension angolo(16,16,1000) open(unit=10,file='distribuzione.dat',status='new') write(6,*)'writing output' do 100 k=1, nd write(6,*)angolo(i,j,k) write(10,*)angolo(i,j,k) 100 continue close(unit=10) end c ------------------ c the subroutine hard is the most important of the program c its contains the dynamics of the hard balls c the input are initial positions (x30,y30,z30) of the first c particle and of the second one (z40,y40,z40), as well as c velocities (v3x0,v3y0,v3z0) for first and second particle c (v4x0,v4y0,v4z0), mass and radius of first (m3,r3) and second c particle (m4,r4) c last but not least, the integer fla specifies the type of c collision: 1 for xe-ag, 2 for xe-ox and 3 for ox-ag c this input is necessary as while oxygen and xenon moves after the c collisio, substrate atoms are not allowed to move, and their mass c and radius parametrize the energy transfer to the substrate and c the corrugation. c Its outputs are final positions and velocities c subroutine hard(x30,y30,z30,x40,y40,z40,v3x0,v3y0,v3z0,v4x0,v4y0,v4z0, % x3,y3,z3,x4,y4,z4,v3x,v3y,v3z,v4x,v4y,v4z,m3,m4,r3,r4,fla) integer fla double precision en3,en2,en30,en20 double precision m3,m4 double precision x3,y3,z3,x4,y4,z4 double precision x30,y30,z30,x40,y40,z40 double precision v3x,v3y,v3z,v4x,v4y,v4z double precision v3x0,v3y0,v3z0,v4x0,v4y0,v4z0 double precision r3,r4 double precision ax,ay,az,bx,by,bz double precision p3x,p3y,p3z,p4x,p4y,p4z double precision p3x0,p3y0,p3z0,p4x0,p4y0,p4z0 double precision t3,t4,t double precision a,b,c,delta,f,check c calculate energies en30=0.5*m3*(v3x0**2+v3y0**2+v3z0**2) en40=0.5*m4*(v4x0**2+v4y0**2+v4z0**2) c calculate momenta p3x0=m3*v3x0 p4x0=m4*v4x0 p3y0=m3*v3y0 p4y0=m4*v4y0 p3z0=m3*v3z0 p4z0=m4*v4z0 c calculate time for collision c relative coordinates a* and velocities b* ax=x30-x40 ay=y30-y40 az=z30-z40 bx=v3x0-v4x0 by=v3y0-v4y0 bz=v3z0-v4z0 c coefficients (a, b, c) of the second order equation whose zeroes c the collision times a=(bx**2+by**2+bz**2) b=2*(ax*bx+ay*by+az*bz) c=ax**2+ay**2+az**2-(r3+r4)**2 c discriminant of the II order equation delta=b**2-4*a*c c cccccccccc if (a.eq.0) then c write(6,*)'3 and 4 have the same velocity: no collision' t=100000000 goto 390 endif c solution of the equation if (delta.lt.0) then c write(6,*)'trajectories do no cross: no collision' t=100000000 c in case of no collision go to 390 and give unaltered velocities c and positions goto 390 else t3=(-b-(delta**0.5))/(2*a) t4=(-b+(delta**0.5))/(2*a) c write(6,*)'t3,t4',t3,t4 c check for errors and identify the first solution in the c future time ( to exclude past collisions) by choosing the less c positive solution. c roundoff errors are excluded c if the lower zero is larger than 1e-17 choose at as time if (t3.gt.(1e-17)) then t=t3 else c otherwise choose the second one t=t4 endif c if also the second one is lower than 1e-17 it means that all c collisions have occurred and the next collision time is at c the infinite future (never) if (t.le.(1e-17)) then t=1000000000 endif endif c write(6,*)'t=',t c ccccccccc c calculate positions at collision 90 x3=x30+v3x0*t y3=y30+v3y0*t z3=z30+v3z0*t x4=x40+v4x0*t y4=y40+v4y0*t z4=z40+v4z0*t check=((x3-x4)**2+(y3-y4)**2+(z3-z4)**2)**0.5 c write(6,*)'check=',check c calculate new momenta c thes eformulas can be demonstrated with some calculations c involving hard balls if (fla.eq.2) then f=(p4x0*(x3-x4)+p4y0*(y3-y4)+p4z0*(z3-z4))/m4- % (p3x0*(x3-x4)+p3y0*(y3-y4)+p3z0*(z3-z4))/m3 f=f/((1/m3+1/m4)*0.5*(r3+r4)**2) p3x=p3x0+f*(x3-x4) p3y=p3y0+f*(y3-y4) p3z=p3z0+f*(z3-z4) p4x=p4x0-f*(x3-x4) p4y=p4y0-f*(y3-y4) p4z=p4z0-f*(z3-z4) c write(6,*)'calculate new momenta collision 1-2' endif c calculate new momenta if ((fla.eq.1).or.(fla.eq.3)) then f=p3x0*(x3-x4)+p3y0*(y3-y4)+p3z0*(z3-z4) f=f/((x3-x4)**2+(y3-y4)**2+(z3-z4)**2) p3x=(x3-x4)*f*(m3-m4)/(m3+m4)+p3x0-(x3-x4)*f p3y=(y3-y4)*f*(m3-m4)/(m3+m4)+p3y0-(y3-y4)*f p3z=(z3-z4)*f*(m3-m4)/(m3+m4)+p3z0-(z3-z4)*f p4x=0 p4y=0 p4z=0 c write(6,*)'calculate new momenta collision with surface' endif v3x=p3x/m3 v3y=p3y/m3 v3z=p3z/m3 v4x=p4x/m4 v4y=p4y/m4 v4z=p4z/m4 en3=0.5*m3*(v3x**2+v3y**2+v3z**2) en4=0.5*m4*(v4x**2+v4y**2+v4z**2) c write(6,*)'initial energies 3,4',en30/1.6e-19,en40/1.6e-19 c write(6,*)'collision position 3',x3/1e-10,y3/1e-10,z3/1e-10 c write(6,*)'collision position 4',x4/1e-10,y4/1e-10,z4/1e-10 c write(6,*)'velocity 3 after',v3x/1e-10,v3y/1e-10,v3z/1e-10 c write(6,*)'velocity 4 after',v4x/1e-10,v4y/1e-10,v4z/1e-10 c write(6,*)'final energies 3,4',en3/1.6e-19,en4/1.6e-19 goto 400 390 x3=x30 y3=y30 z3=z30 v3x=v3x0 v3y=v3y0 v3z=v3z0 x4=x40 y4=y40 z4=z40 v4x=v4x0 v4y=v4y0 v4z=v4z0 400 end c ------------------- c the subroutine tim(...) includes the kinematics. i.t. it calculate c from initial positions of first (x30,y30,z30) an second c partcile (x40,y40,z40) and velocities (v3x0,v3y0,v3z0) and c (v4x0,v4y0,v4z0), the time t of the next collision between them c As the particle have a size it requries as input also their radii c the time of impact is determined as the time at whch the distance c of their centers of mass equals the sum of their radii c (particles are always assumed to be spheres) subroutine tim(x30,y30,z30,x40,y40,z40,v3x0,v3y0,v3z0,v4x0,v4y0,v4z0, % r3,r4,t) double precision en3,en2,en30,en20 double precision x30,y30,z30,x40,y40,z40 double precision v3x0,v3y0,v3z0,v4x0,v4y0,v4z0 double precision r3,r4 double precision ax,ay,az,bx,by,bz double precision t3,t4,t double precision a,b,c,delta,f,check c write(6,*)'in sub time',x30,y30,z30 c write(6,*)'in sub time',v3x0,v3y0,v3z0 c calculate tiem for colliison ax=x30-x40 ay=y30-y40 az=z30-z40 bx=v3x0-v4x0 by=v3y0-v4y0 bz=v3z0-v4z0 c coefficients of second order equation whose zeros are collision c times a=(bx**2+by**2+bz**2) b=2*(ax*bx+ay*by+az*bz) c=ax**2+ay**2+az**2-(r3+r4)**2 c determinant of the equation delta=b**2-4*a*c if (a.eq.0) then c write(6,*)'3 and 4 have the same velocity: no collision' t=100000000 goto 390 endif if (delta.lt.0) then c write(6,*)'trajectories do no cross: no collision' t=100000000 goto 390 else t3=(-b-(delta**0.5))/(2*a) t4=(-b+(delta**0.5))/(2*a) c write(6,*)'t3,t4',t3,t4 if (t3.gt.(1e-17)) then t=t3 else t=t4 endif if (t.le.(1e-17)) then t=1000000000 endif endif c write(6,*)'t=',t 390 end --------------- Here is a typical input file (simuloluca_bur.dat) --------------- 0.5 5.5 4 0 60 3 1000 250 0.3 0.1 0.4 prova_bur2.dat 131 32 2.2 1.125 2.88 100 dprova_bur2.dat ------ the output file (prova_bur2.dat) looks as follows: ------------------------------- program sigmaluca_bur_rev.for hard cube model During the simulation some particles may be lost as they go out of the lattice if such effect is larger than 2 % an ALERT appears m1,m2,m (amu)= 131.000000000000 32.0000000000000 250.000000000000 r1,r2,r,sep=(ang) 2.20000000000000 1.12500000000000 3.60600000000000 0.100000000000000 corr= 0.300000000000000 zox=(ang) 0.774794578835283 nr= 16 npart= 1000 JJ LL ENERGY TETA SIGMA DSIGMA SIGMA-DIR DSIGMA-DIR avexe(jj,ll),aveox(jj,ll),avesurf(jj,ll) avexedes(jj,ll),aveoxdes(jj,ll),avesurfdes(jj,ll) avexedesn(jj,ll),aveoxdesn(jj,ll) 1 1 0.500000000000000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 7.111532156512902E-002 7.166089971251584E-003 0.421718596576044 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 ERROR on total particles numbers (defect) ndes,nrim,oxins,TOT 0 998 0 998 particle no colliding= 0 particle lost 0 1 2 0.500000000000000 30.0000000000000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.166930158228453 1.053821225784938E-002 0.322531657684004 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 ERROR on total particles numbers (defect) ndes,nrim,oxins,TOT 0 998 0 998 particle no colliding= 0 particle lost 0 1 3 0.500000000000000 60.0000000000000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.345668217100521 2.094979037909959E-002 0.133382020690687 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 ERROR on total particles numbers (defect) ndes,nrim,oxins,TOT 1 996 0 997 particle no colliding= 0 particle lost 0 2 1 2.16666666666667 0.000000000000000E+000 0.401203602552414 0.401203602552414 0.401203602552414 0.401203602552414 0.312372359821025 3.211911304327007E-002 1.82217526958237 0.456720099579742 0.854569107222673 0.855377459864252 6.779504727735482E-002 0.415719099586430 ERROR on total particles numbers (defect) ndes,nrim,oxins,TOT 13 985 0 998 particle no colliding= 0 particle lost 0 2 2 2.16666666666667 30.0000000000000 5.21042108535767 1.44511067867279 5.21042108535767 1.44511067867279 0.728717005319525 5.180487215783579E-002 1.38614491126066 0.273863713904288 0.815530662590651 1.07727229017173 8.675670645132080E-002 0.485753677381981 ERROR on total particles numbers (defect) ndes,nrim,oxins,TOT 14 983 0 997 particle no colliding= 0 particle lost 0 2 3 2.16666666666667 60.0000000000000 5.61685037612915 1.50116658210754 4.81444358825684 1.38981008529663 1.50769228570874 8.579652677242114E-002 0.573177929965491 0.596901395428841 1.04955176496075 0.520213506277073 5.808792640565905E-002 0.503620404069115 ERROR on total particles numbers (defect) ndes,nrim,oxins,TOT 25 974 0 999 particle no colliding= 0 particle lost 0 3 1 3.83333333333333 0.000000000000000E+000 10.0100097656250 2.00200200080872 10.0100097656250 2.00200200080872 0.563364257059287 8.604808087963160E-002 3.18392097569962 0.894564278965764 1.30056950334514 1.63819955102243 0.180228528972153 0.659641658352532 ERROR on total particles numbers (defect) ndes,nrim,oxins,TOT 21 978 0 999 particle no colliding= 0 particle lost 0 3 2 3.83333333333333 30.0000000000000 8.40840816497803 1.83486509323120 8.00800800323486 1.79064500331879 1.26347318455839 0.102748281280159 2.46711184779998 0.593155908113280 1.73511179641651 1.50506562880354 8.061656016461281E-002 0.763099954418091 3 3 3.83333333333333 60.0000000000000 16.3999996185303 2.56124973297119 14.0000000000000 2.36643195152283 2.66895984319189 0.147436119403685 1.01693743293297 1.08055483652816 1.86915372265276 0.883624774152415 0.106083151018660 0.754689979528098 ERROR on total particles numbers (defect) ndes,nrim,oxins,TOT 28 971 0 999 particle no colliding= 0 particle lost 0 4 1 5.50000000000000 0.000000000000000E+000 11.2112112045288 2.11871981620789 11.2112112045288 2.11871981620789 0.799751397059480 9.580639652239514E-002 4.60444217816035 1.46769020722464 1.56709138202956 2.46521841074580 0.235125270936545 0.783773105254253 ERROR on total particles numbers (defect) ndes,nrim,oxins,TOT 33 965 0 998 particle no colliding= 0 particle lost 0 4 2 5.50000000000000 30.0000000000000 13.2264528274536 2.30242991447449 13.2264528274536 2.30242991447449 1.80207097209788 0.126831423788442 3.57109791398705 1.23012372789872 2.08067027985833 2.18920599224296 0.249712342749198 0.869781999511122 ERROR on total particles numbers (defect) ndes,nrim,oxins,TOT 56 942 0 998 particle no colliding= 0 particle lost 0 4 3 5.50000000000000 60.0000000000000 22.4448909759521 2.99932479858398 18.0360717773438 2.68865871429443 3.80873563599367 0.206748792980248 1.48451588089947 1.82121687734343 2.46433569517776 1.21444742747881 0.225001848709381 0.930735299621640 -------------- **************************************************************************************************** //////////////////////////////////////////////////////////////////////////////////////////////////// **************************************************************************************************** Listing 9 - FORTRAN short: program sigmaluca_bur_rev double precision tetamin,tetamax double precision x10,y10,z10 double precision x1,y1,z1 double precision x20,y20,z20 double precision x2,y2,z2 double precision x3,y3,z3,x4,y4, % z4,v3x, % v3y,v3z,v4x,v4y,v4z double precision v1x0,v1y0,v1z0 double precision v1x,v1y,v1z double precision v2x0,v2y0,v2z0 double precision v2x,v2y,v2z double precision x1in,y1in,z1in,dt, % en,en2,teta double precision x0,y0,z0 double precision x,y,z double precision diffen,enin,enfin, % en10,en20,en1norm,en2norm double precision avexedes,avexedesn, % aveoxdes,aveoxdesn double precision tetaout character*20 luca1,luca2 integer ndes,nrim,nx,ny,xeins,oxins, % flad,ndir,nteta,index integer nocoll,nlost double precision sigmadir,deltasigmadir double precision vx0,vy0,vz0 double precision vx,vy,vz double precision xdes,ydes,deltasigmades double precision emin,emax,sigmades,ubin, % x2in,y2in,z2in,zmin integer nr,flag,i1ag,j1ag,i2ag,j2ag, % npart,nen double precision aveox,avexe, % avesurf,avesurfdes double precision t,t12,t1ag,t2ag,tempo double precision en1x,en1y,en1z double precision r1,r2,r,a,sep double precision m1,m2,m,corr dimension aveox(20,20),avexe(20,20), % avesurf(20,20) dimension avexedes(20,20),aveoxdes(20,20), % avesurfdes(20,20) dimension avexedesn(20,20), % aveoxdesn(20,20) dimension en(20,20),teta(20,20) dimension sigmades(20,20),sigmadir(20,20) dimension deltasigmades(20,20), % deltasigmadir(20,20) dimension xdes(1000),ydes(1000) dimension x(16,16),y(16,16),z(16,16) dimension x0(16,16),y0(16,16),z0(16,16) dimension vx0(16,16),vy0(16,16),vz0(16,16) dimension vx(16,16),vy(16,16),vz(16,16) dimension tetaout(16,16,1000) common/dati/ x10,y10,z10,x20,y20, % z20,x0,y0,z0,v1x0,v1y0,v1z0, % v2x0,v2y0,v2z0,vx0,vy0,vz0,m1,m2,m,r1,r2,r nx=16523 ny=98223 open(unit=2,file='simu.dat',status='old') read(2,*)emin,emax,nen read(2,*)tetamin,tetamax,nteta read(2,*)npart read(2,*)m,corr,sep read(2,*)ubin read(2,5)luca1 read(2,*)m1,m2 read(2,*)r1,r2 read(2,*)a read(2,*)index read(2,5)luca2 5 format(A) close(unit=2) r=0.5*((a/2)**2/corr+corr) m1=m1*1.67e-27 m2=m2*1.67e-27 m=m*1.67e-27 r=r*1e-10 sep=sep*1e-10s r1=r1*1e-10 r2=r2*1e-10 nr=16 a=a*1e-10 x2in=0 y2in=0 if (index.eq.100) then z2in=((r+r2+sep)**2- & (a/2*(2**0.5))**2)**0.5-r endif if (index.eq.111) then z2in=((r+r2+sep)**2- % (a/4*(3**0.5))**2)**0.5-r endif if (index.eq.1001) then z2in=((r+r2+sep)**2-(a/2)**2)**0.5-r endif if (index.eq.1002) then z2in=((r+r2+sep)**2-(a/2)**2)**0.5-r endif open(unit=3,file=luca1,status='new') write(3,*)'m1,m2,m (amu)=',m1/1.67e-27, % m2/1.67e-27,m/1.67e-27 write(3,*)'r1,r2,r,sep=(ang)',r1/1e-10, % r2/1e-10,r/1e-10,sep/1e-10 write(3,*)'corr=',corr write(3,*)'zox=(ang)',z2in/1e-10 write(3,*)'nr=',nr write(3,*)'npart=',npart write(3,*)'JJ LL ENERGY TETA' write(3,*)' SIGMA DSIGMA' write(3,*)' SIGMA-DIR DSIGMA-DIR' write(3,*)'avexe(jj,ll),aveox(jj,ll), % avesurf(jj,ll)' write(3,*)'avexedes(jj,ll),aveoxdes(jj,ll), % avesurfdes(jj,ll)' write(3,*)'avexedesn(jj,ll),aveoxdesn(jj,ll)' write(6,*)'m1,m2,m=',m1/1.67e-27, % m2/1.67e-27,m/1.67e-27 write(6,*)'r1,r2,r,sep=',r1/1e-10, % r2/1e-10,r/1e-10,sep/1e-10 write(6,*)'corr=',corr write(6,*)'zox=(ang)',z2in/1e-10 write(6,*)'nr=',nr write(6,*)'npart=',npart write(6,*)2 open(unit=10,file=luca2,status='new') do 20000 jj=1, nen do 15000 ll=1, nteta if (nen.eq.1) then en(1,ll)=emin else en(jj,ll)=emin+(emax-emin)*(jj-1)/(nen-1) endif if (nteta.eq.1) then teta(jj,1)=tetamin else teta(jj,ll)=tetamin+(tetamax- % tetamin)*(ll-1)/(nteta-1) endif ndes=0 ndir=0 nrim=0 nocoll=0 nlost=0 xeins=0 oxins=0 sigmades(jj,ll)=0 sigmadir(jj,ll)=0 aveox(jj,ll)=0 avesurfdes(jj,ll)=0 aveoxdes(jj,ll)=0 avexedes(jj,ll)=0 aveoxdesn(jj,ll)=0 avexedesn(jj,ll)=0 avexe(jj,ll)=0 avesurf(jj,ll)=0 do 10000 ii=1, npartcle if ((ndes+nrim).gt.(npart)) then write(3,*)'ERROR on total particles numbers' goto 10010 endif en1y=0 en1x=en(jj,ll)*(sin(3.1415* % teta(jj,ll)/180.))**2 en1z=en(jj,ll)*(cos(3.1415* % teta(jj,ll)/180.))**2 z10=100 x10=-10+20*ran(nx) y10=-10+20*ran(ny) x1in=x10 y1in=y10 x10=x10-z10*tan(3.1415*teta(jj,ll)/180) x10=x10*1e-10 y10=y10*1e-10 z10=z10*1e-10 z1in=z10 x1=x10 y1=y10 z1=z10 x20=x2in y20=y2in z20=z2in x2=x20 y2=y20 z2=z20 do 10 i=1, nr do 9 j=1, nr if (index.eq.100) then x0(i,j)=-0.5*(nr-1)*a+(i-1)*a y0(i,j)=-0.5*(nr-1)*a+(j-1)*a z0(i,j)=-r endif if (index.eq.111) then x0(i,j)=-0.5*a+(i-1)*a+(j-1)*a/2- % (nr*0.5-1)*a-(nr*0.5-1)*a/2 y0(i,j)=a*(3.**0.5)/4+(j-1)*a*(3.**0.5)/2- % (nr*0.5-1)*a*(3.**0.5)/2 z0(i,j)=-r endif if (index.eq.1001) then x0(i,j)=-0.5*(nr-1)*a+(i-1)*a y0(i,j)=-0.5*(nr-1)*a+(j-1)*a+0.5*a z0(i,j)=-r endif if (index.eq.1002) then x0(i,j)=-0.5*(nr-1)*a+(i-1)*a+0.5*a y0(i,j)=-0.5*(nr-1)*a+(j-1)*a z0(i,j)=-r endif x(i,j)=x0(i,j) y(i,j)=y0(i,j) z(i,j)=z0(i,j) 9 continue 10 continue en1x=en1x*1.6e-19 en1y=en1y*1.6e-19 en1z=en1z*1.6e-19 v1x0=(2*en1x/m1)**0.5 v1y0=(2*en1y/m1)**0.5 v1z0=-(2*en1z/m1)**0.5 v1x=v1x0 v1y=v1y0 v1z=v1z0 v2x0=0 v2y0=0 v2z0=0 v2x=v2x0 v2y=v2y0 v2z=v2z0 do 20 i=1, nr do 19 j=1, nr vx0(i,j)=0 vy0(i,j)=0 vz0(i,j)=0 vx(i,j)=0 vy(i,j)=0 vz(i,j)=0 19 continue 20 continue tempo=0 en20=0.5*m2*(v2x0**2+v2y0**2+v2z0**2) en10=0.5*m1*(v1x0**2+v1y0**2+v1z0**2) tempo=0 flad=0 do 1000 kk=1, 300 t1ag=1000000 flag=0 dt=0 do 100 i=1, nr do 99 j=1, nr call tim(x10,y10,z10,x0(i,j), % y0(i,j),z0(i,j),v1x0,v1y0,v1z0, % vx0(i,j),vy0(i,j),vz0(i,j),r1,r,t) if (t.lt.t1ag) then t1ag=t i1ag=i j1ag=j endif 99 continue 100 continue t12=1000000 call tim(x10,y10,z10,x20,y20,z20, % v1x0,v1y0,v1z0,v2x0,v2y0,v2z0, % r1,r2,t) if (t.lt.t12) then t12=t endif t2ag=1000000 do 200 i=1, nr do 199 j=1, nr call tim(x20,y20,z20,x0(i,j), % y0(i,j),z0(i,j),v2x0,v2y0,v2z0, % vx0(i,j),vy0(i,j),vz0(i,j),r2,r,t) if (t.lt.t2ag) then t2ag=t i2ag=i j2ag=j endif 199 continue 200 continue if ((t1ag.lt.t2ag).and.(t1ag.lt.t12)) then flag=1 dt=t1ag tempo=tempo+dt endif if ((t12.lt.t1ag).and.(t12.lt.t2ag)) then flag=2 dt=t12 tempo=tempo+dt endif if ((t2ag.lt.t1ag).and.(t2ag.lt.t12)) then flag=3 dt=t2ag tempo=tempo+dt endif if (kk.eq.1) then flad=flag endif if ((dt.gt.10000).and.(kk.eq.1)) then flad=0 endif if (flag.eq.0) then dt=1e-10 kk=299 endif if (flad.eq.0) then nocoll=nocoll+1 endif x1=x10+v1x0*dt y1=y10+v1y0*dt z1=z10+v1z0*dt x2=x20+v2x0*dt y2=y20+v2y0*dt z2=z20+v2z0*dt if (flag.eq.0) then v1x=v1x0 v1y=v1y0 v1z=v1z0 v2x=v2x0 v2y=v2y0 v2z=v2z0 goto 290 endif if (flag.eq.1) then call hard(x10,y10,z10,x0(i1ag,j1ag) % ,y0(i1ag,j1ag),z0(i1ag,j1ag), % v1x0,v1y0,v1z0,vx0(i1ag,j1ag), % vy0(i1ag,j1ag),vz0(i1ag,j1ag), % x3,y3,z3,x4,y4,z4,v3x,v3y,v3z,v4x, % v4y,v4z,m1,m,r1,r,flag) v1x=v3x v1y=v3y v1z=v3z v2x=v2x0 v2y=v2y0 v2z=v2z0 goto 290 endif if (flag.eq.2) then call hard(x10,y10,z10,x20,y20,z20, % v1x0,v1y0,v1z0,v2x0,v2y0,v2z0, % x3,y3,z3,x4,y4,z4,v3x,v3y,v3z,v4x, % v4y,v4z,m1,m2,r1,r2,flag) v1x=v3x v1y=v3y v1z=v3z v2x=v4x v2y=v4y v2z=v4z goto 290 endif if (flag.eq.3) then call hard(x20,y20,z20,x0(i2ag,j2ag), % y0(i2ag,j2ag),z0(i2ag,j2ag), % v2x0,v2y0,v2z0,vx0(i2ag,j2ag), % vy0(i2ag,j2ag),vz0(i2ag,j2ag), % x3,y3,z3,x4,y4,z4,v3x,v3y,v3z, % v4x,v4y,v4z,m2,m,r2,r,flag) v1x=v1x0 v1y=v1y0 v1z=v1z0 v2x=v3x v2y=v3y v2z=v3z goto 290 endif 290 if (z1.lt.0) then xeins=xeins+1 goto 10000 endif if (z2.lt.0) then oxins=oxins+1 goto 10000 endif x10=x1 y10=y1 z10=z1 x20=x2 y20=y2 z20=z2 v1x0=v1x v1y0=v1y v1z0=v1z v2x0=v2x v2y0=v2y v2z0=v2z v10=(v1x0**2+v1y0**2+v1z0**2)**0.5 v20=(v2x0**2+v2y0**2+v2z0**2)**0.5 if (v10.eq.0) then goto 900 else endif 900 if (v20.eq.0) then goto 902 else endif 902 if (z10.gt.z1in) then goto 1001 endif if (flad.eq.0) then nlost=nlost+1 write(3,*)'LOST PARTICLE' write(3,*)'1 coordinate ',x10,y10,z10 write(3,*)'1 velocity',v1x0,v1y0,v1z0 write(3,*)'2 coordinate ',x20,y20,z20 write(3,*)'2 velocity',v2x0,v2y0,v2z0 write(3,*)'t12,t1ag,t2ag,t',t12,t1ag,t2ag,t goto 1001 endif 1000 continue 1001 en2norm=0.5*m2*v2z**2 en1norm=0.5*m1*v1z**2 en1=0.5*m1*(v1x**2+v1y**2+v1z**2) en2=0.5*m2*(v2x**2+v2y**2+v2z**2) aveox(jj,ll)=aveox(jj,ll)+en2 avexe(jj,ll)=avexe(jj,ll)+en1 avesurf(jj,ll)=avesurf(jj,ll) % +(en10-en1-en2) if ((en2norm.gt.(ubin*1.6e-19)). % and.(v2z.gt.0)) then ndes=ndes+1 if (flad.eq.2) then ndir=ndir+1 endif write(6,*)'ox velocit dir', % v2x0/v20,v2y0/v20,v2z0/v20 tetaout(jj,ll,ndes)=180./3.1415* % acos(v2z/(v2x**2+v2y**2+v2z**2)**0.5) write(6,*)'tetaout=',v2z/ % (v2x**2+v2y**2+v2z**2)**0.5, % tetaout(jj,ll,ndes) avesurfdes(jj,ll)=avesurfdes(jj,ll) % +(en(jj,ll)*1.6e-19)-0.5*m2 % *(v2x**2+v2y**2+v2z**2)- % 0.5*m1*(v1x**2+v1y**2+v1z**2) aveoxdes(jj,ll)=aveoxdes(jj,ll)+ % 0.5*m2*(v2x**2+v2y**2+v2z**2) aveoxdesn(jj,ll)=aveoxdesn(jj,ll) % +0.5*m2*(v2z**2) avexedes(jj,ll)=avexedes(jj,ll) % +0.5*m1*(v1x**2+v1y**2+v1z**2) avexedesn(jj,ll)=avexedesn(jj,ll) % +0.5*m1*(v1z**2) xdes(ndes)=x1in ydes(ndes)=y1in else nrim=nrim+1 endif 10000 continue if ((ndes+nrim).gt.(npart)) then write(6,*)'ERROR on total particles numbers (excess)' write(6,*)'ndes,nrim,ndes+nrim',ndes,nrim,ndes+nrim write(3,*)'ERROR on total particles numbers (excess)' write(3,*)'ndes,nrim,ndes+nrim',ndes,nrim,ndes+nrim write(10,*)'ERROR on total particles numbers (excess)' write(10,*)'ndes,nrim,ndes+nrim',ndes,nrim,ndes+nrim endif if ((ndes+nrim+oxins+nocoll+nlost).lt.(npart)) then write(6,*)'ERROR on total particles numbers (defect)' write(6,*)'ndes,nrim,oxins,TOT',ndes,nrim,oxins,ndes+nrim+oxins write(3,*)'ERROR on total particles numbers (defect)' write(3,*)'ndes,nrim,oxins,TOT',ndes,nrim,oxins,ndes+nrim+oxins write(10,*)'ERROR on total particles numbers (defect)' write(10,*)'ndes,nrim,oxins,TOT',ndes,nrim,oxins,ndes+nrim+oxins write(6,*)'particle no colliding=',nocoll write(3,*)'particle no colliding=',nocoll write(10,*)'particle no colliding=',nocoll write(6,*)'particle lost',nlost write(3,*)'particle lost',nlost write(10,*)'particle lost',nlos endif 10010 sigmades(jj,ll)=ndes*1.0/((ndes+nrim)*1.0)*400.0 sigmadir(jj,ll)=ndir*1.0/((ndes+nrim)*1.0)*400.0 if (((nocoll+oxins+xeins+nlost)*1.0).gt.(0.02*npart)) then write(6,*)'nocoll=',nocoll write(6,*)'oxins=',oxins write(6,*)'xeins=',xeins write(6,*)'ndes,nrim=',ndes,nrim write(6,*)'npart=',npart write(6,*)' ALERT ! TOO PARTICLES INSIDE OR LOST' write(6,*)' TAKE RESULT WITH CARE POSSIBLE FLAW !!' write(3,*)'nocoll=',nocoll write(3,*)'oxins=',oxins write(3,*)'xeins=',xeins write(3,*)'ndes,nrim=',ndes,nrim write(3,*)'nlost=',nlost write(3,*)'npart=',npart write(3,*)(nocoll+oxins+xeins+nlost)*1.0 write(3,*)(0.02*npart) write(3,*)' ALERT ! TOO PARTICLES INSIDE OR LOST' write(3,*)' TAKE RESULT WITH CARE POSSIBLE FLAW !!' write(10,*)'nocoll=',nocoll write(10,*)'oxins=',oxins write(10,*)'xeins=',xeins write(10,*)'ndes,nrim=',ndes,nrim write(10,*)'npart=',npart write(10,*)' ALERT ! TOO PARTICLES INSIDE OR LOST' write(10,*)' TAKE RESULT WITH CARE POSSIBLE FLAW !!' write(6,*)'Penetration in the bulk may indicate that there' write(6,*)'are holes between the atoms which model the surface' write(6,*)'Lost particles may indicate that the set of initial' write(6,*)'parameters is too wide for the dimension of the' write(6,*)'array of atoms so that some particle never reach it' write(3,*)'Penetration in the bulk may indicate that there' write(3,*)'are holes between the atoms which model the surface' write(3,*)'Lost particles may indicate that the set of initial' write(3,*)'parameters is too wide for the dimension of the' write(3,*)'array of atoms so that some particle never reach it' endif deltasigmades(jj,ll)=(ndes**0.5)*1.0/((ndes+nrim)*1.0)*400.0 deltasigmadir(jj,ll)=(ndir**0.5)*1.0/((ndes+nrim)*1.0)*400.0 aveox(jj,ll)=aveox(jj,ll)/((ndes+nrim)*1.6e-19) avexe(jj,ll)=avexe(jj,ll)/((ndes+nrim)*1.6e-19) avesurf(jj,ll)=avesurf(jj,ll)/((ndes+nrim)*1.6e-19) if (ndes.gt.0) then aveoxdes(jj,ll)=aveoxdes(jj,ll)/(ndes*1.0)/1.6e-19 aveoxdesn(jj,ll)=aveoxdesn(jj,ll)/(ndes*1.0)/1.6e-19 avexedes(jj,ll)=avexedes(jj,ll)/(ndes*1.0)/1.6e-19 avexedesn(jj,ll)=avexedesn(jj,ll)/(ndes*1.0)/1.6e-19 avesurfdes(jj,ll)=avesurfdes(jj,ll)/(ndes*1.0)/1.6e-19 else aveoxdes(jj,ll)=0 aveoxdesn(jj,ll)=0 avexedes(jj,ll)=0 avexedesn(jj,ll)=0 avesurfdes(jj,ll)=0 endif write(10,*)jj,ll,en(jj,ll),teta(jj,ll),ndes if (ndes.gt.0) then call conto(jj,ll,ndes,tetaout) else write(10,*) endif write(10,*)'----------' write(6,*)jj,ll,en(jj,ll),teta(jj,ll) write(6,*)sigmades(jj,ll),deltasigmades(jj,ll) write(6,*)sigmadir(jj,ll),deltasigmadir(jj,ll) write(6,*)'---------------' write(3,*)jj,ll,en(jj,ll),teta(jj,ll) write(3,*)sigmades(jj,ll),deltasigmades(jj,ll) write(3,*)sigmadir(jj,ll),deltasigmadir(jj,ll) write(3,*)avexe(jj,ll),aveox(jj,ll),avesurf(jj,ll) write(3,*)avexedes(jj,ll),aveoxdes(jj,ll),avesurfdes(jj,ll) write(3,*)avexedesn(jj,ll),aveoxdesn(jj,ll) write(3,*) 15000 continue 20000 continue close(unit=3) close(unit=10) stop end subroutine conto(i,j,n,angolo) integer i,j,n,nt double precision angolo,t,tmin,tmax dimension angolo(16,16,1000),nt(9) do 10 l=1, 9 nt(l)=0 10 continue do 100 k=1, n do 80 l=1, 9 tmin=(l-1)*10. tmax=l*10. if ((angolo(i,j,k).gt.tmin).and. % (angolo(i,j,k).lt.tmax)) then nt(l)=nt(l)+1 endif 80 continue 100 continue do 110 l=1, 9 write(10,*)(l-0.5)*10.,nt(l)/(n*1.0) 110 continue end subroutine hard(x30,y30,z30,x40, % y40,z40,v3x0,v3y0,v3z0,v4x0,v4y0, % v4z0,x3,y3,z3,x4,y4,z4,v3x,v3y, % v3z,v4x,v4y,v4z,m3,m4,r3,r4,fla) integer fla double precision en3,en2,en30,en20 double precision m3,m4 double precision x3,y3,z3,x4,y4,z4 double precision x30,y30,z30,x40, % y40,z40 double precision v3x,v3y,v3z,v4x, % v4y,v4z double precision v3x0,v3y0,v3z0, % v4x0,v4y0,v4z0 double precision r3,r4 double precision ax,ay,az,bx,by,bz double precision p3x,p3y,p3z,p4x, % p4y,p4z double precision p3x0,p3y0,p3z0, % p4x0,p4y0,p4z0 double precision t3,t4,t double precision a,b,c,delta,f,check en30=0.5*m3*(v3x0**2+v3y0**2+v3z0**2) en40=0.5*m4*(v4x0**2+v4y0**2+v4z0**2) p3x0=m3*v3x0 p4x0=m4*v4x0 p3y0=m3*v3y0 p4y0=m4*v4y0 p3z0=m3*v3z0 p4z0=m4*v4z0 ax=x30-x40 ay=y30-y40 az=z30-z40 bx=v3x0-v4x0 by=v3y0-v4y0 bz=v3z0-v4z0 a=(bx**2+by**2+bz**2) b=2*(ax*bx+ay*by+az*bz) c=ax**2+ay**2+az**2-(r3+r4)**2 delta=b**2-4*a*c if (a.eq.0) then t=100000000 goto 390 endif if (delta.lt.0) then t=100000000 goto 390 else t3=(-b-(delta**0.5))/(2*a) t4=(-b+(delta**0.5))/(2*a) if (t3.gt.(1e-17)) then t=t3 else t=t4 endif if (t.le.(1e-17)) then t=1000000000 endif endif 90 x3=x30+v3x0*t y3=y30+v3y0*t z3=z30+v3z0*t x4=x40+v4x0*t y4=y40+v4y0*t z4=z40+v4z0*t check=((x3-x4)**2+(y3-y4)**2+ % (z3-z4)**2)**0.5 if (fla.eq.2) then f=(p4x0*(x3-x4)+p4y0*(y3-y4)+ % p4z0*(z3-z4))/m4- % (p3x0*(x3-x4)+p3y0*(y3-y4)+ % p3z0*(z3-z4))/m3 f=f/((1/m3+1/m4)*0.5*(r3+r4)**2) p3x=p3x0+f*(x3-x4) p3y=p3y0+f*(y3-y4) p3z=p3z0+f*(z3-z4) p4x=p4x0-f*(x3-x4) p4y=p4y0-f*(y3-y4) p4z=p4z0-f*(z3-z4) endif if ((fla.eq.1).or.(fla.eq.3)) then f=p3x0*(x3-x4)+p3y0*(y3-y4)+ % p3z0*(z3-z4) f=f/((x3-x4)**2+(y3-y4)**2+ % (z3-z4)**2) p3x=(x3-x4)*f*(m3-m4)/(m3+m4)+ % p3x0-(x3-x4)*f p3y=(y3-y4)*f*(m3-m4)/(m3+m4)+ % p3y0-(y3-y4)*f p3z=(z3-z4)*f*(m3-m4)/(m3+m4)+ % p3z0-(z3-z4)*f p4x=0 p4y=0 p4z=0 endif v3x=p3x/m3 v3y=p3y/m3 v3z=p3z/m3 v4x=p4x/m4 v4y=p4y/m4 v4z=p4z/m4 en3=0.5*m3*(v3x**2+v3y**2+v3z**2) en4=0.5*m4*(v4x**2+v4y**2+v4z**2) goto 400 390 x3=x30 y3=y30 z3=z30 v3x=v3x0 v3y=v3y0 v3z=v3z0 x4=x40 y4=y40 z4=z40 v4x=v4x0 v4y=v4y0 v4z=v4z0 400 end subroutine tim(x30,y30,z30,x40,y40, % z40,v3x0,v3y0,v3z0,v4x0,v4y0,v4z0, % r3,r4,t) double precision en3,en2,en30,en20 double precision x30,y30,z30,x40, % y40,z40 double precision v3x0,v3y0,v3z0, % v4x0,v4y0,v4z0 double precision r3,r4 double precision ax,ay,az,bx,by,bz double precision t3,t4,t double precision a,b,c,delta,f,check ax=x30-x40 ay=y30-y40 az=z30-z40 bx=v3x0-v4x0 by=v3y0-v4y0 bz=v3z0-v4z0 a=(bx**2+by**2+bz**2) b=2*(ax*bx+ay*by+az*bz) c=ax**2+ay**2+az**2-(r3+r4)**2 delta=b**2-4*a*c if (a.eq.0) then t=100000000 goto 390 endif if (delta.lt.0) then t=100000000 goto 390 else t3=(-b-(delta**0.5))/(2*a) t4=(-b+(delta**0.5))/(2*a) if (t3.gt.(1e-17)) then t=t3 else t=t4 endif if (t.le.(1e-17)) then t=1000000000 endif endif 390 end **************************************************************************************************** //////////////////////////////////////////////////////////////////////////////////////////////////// **************************************************************************************************** Listing A5: #include // standard "C" header #include // for simple mathematical functions #include // for simple input and output operations #include "nr.h" // from numerical recipes #define NRANSI // some compiler require this statement, // which indicates that standard "C" code follows int idum=(-1); // for the random number generator, include always this line in the examples // KMCS, modified Kisliuk model, conventional programming ste/bgh //**** model parameters ********************************************* int nx = 50; // size of the lattice without boundary (square lattice) int Average = 5; // average results of the last "average" Monte Carlo steps float SaturationCoverage = 0.8; // iteration will be stopped if satcov is reached float Flux= 1.0; // flux of the impinging particles number of particles per sec and normalized to // the size of the grid i.e. Flux = 2500 on a 50 by 50 grid equals Flux = 1.0 ML/s float Pclean = 0.6; // 1 - Sclean, desorption probability from clean site float Poccu = 0.7; // 1 - Soccu, desorption probability from occupied site float Pdiff = 1.2; // diffusion probability /*** program variables, do not change ********************************/ int ny; float Sclean; // adsorption probability on clean site float Soccu; // adsorption probability on occupied site float AdsorptionProbability[200]; // calculated adsorption probability float SAverage[200]; // increase indexing if necessary float Coverage[200]; // calculated coverage float RealTime[200]; // calculated exposure time int Lattice[200][200]; // =1 lattice site occupied, =0 not occupied int Precursor[200][200]; // =1 precursor state occupied, =0 not occupied // the linear selection procedure requires to calculate a difference of numbers which can have // very similar numerical values. therefore we use double precision for those variables double Desorption[200][200]; // every site on the surface will be assigned to double Deposition[200][200]; // a desorption, deposition, and diffusion probability double Diffusion[200][200]; int Neighbors[200][200]; // number of occupied nearest neighbor of site i,j double CurrentTime; double OldTime; int Sum; // number of adsorbed particles double Omega; // total probability of all processes for the given configuration double Depos; // total probability of deposition for the given configuration double dt; // time increment double K; // used for the KMCS selection procedure int index; // for storing the data char answer[1]; void CalculateOmega(int i, int j) { // calculates omega and assigns all probabilities according to the precursor rules { Omega = Omega - (Desorption[i][j]+Diffusion[i][j]); Depos = Depos - Deposition[i][j]; Neighbors[i][j]=0; if (Precursor[i][j-1]==1) Neighbors[i][j] ++; if (Precursor[i][j+1]==1) Neighbors[i][j] ++; if (Precursor[i+1][j]==1) Neighbors[i][j] ++; if (Precursor[i-1][j]==1) Neighbors[i][j] ++; if (Precursor[i][j]==0) { // clean precursor site... Desorption[i][j]=0.0; Diffusion[i][j]=0.0; if (Lattice[i][j]==0) Deposition[i][j]=Sclean; // and clean surface site else Deposition[i][j]=Soccu; // and occupied surface site } else { // precursor occupied Desorption[i][j]=Poccu; Deposition[i][j]=0.0; Diffusion[i][j]=(4.0-Neighbors[i][j])*Pdiff; } Omega = Omega + (Desorption[i][j]+Diffusion[i][j]); Depos = Depos + Deposition[i][j]; } void Initialize() { int i,j; // loop index, adsorption site Omega=0; // total transition probability of the configuration Depos=0; // probability for deposition/adsorption of the configuration for (i=1; i<=nx; i++) { for (j=1; j<=ny; j++) { Lattice[i][j]=0; Precursor[i][j]=0; } } for (i=0; i<=99; i++) { AdsorptionProbability[i]=0; Coverage[i]=0; SAverage[i]=0; RealTime[i]=0; } for (i=1; i<=nx; i++) { for (j=1; j<=ny; j++) { Desorption[i][j] = 0.0; Deposition[i][j] = 0.0; Diffusion[i][j] = 0.0; CalculateOmega(i,j); } } } void DepositionProcess() // deposition i.e. adsorption of particles { int Flag; // = 0 deposition successful, = 1 no deposition int i,j; // adsorption site Flag = 1; if ( K < (Depos*Flux/(nx*ny)) ) { for (i=2; ((i<=nx-2) & (Flag == 1)); i++) { // go through the lattice until an adsorption took place for (j=2; ((j<=ny-2) & (Flag == 1)); j++) { // go through the lattice until an adsorption took place if (Flag==1) K = K - Deposition[i][j]*Flux/(nx*ny); // linear selection rule if ((K <=0) & (Flag==1)) { // deposition has been chosen as the process to proceed Flag = 0; if (Lattice[i][j]==0) { // adsorption on free site Lattice[i][j]=1; Sum++; CalculateOmega(i, j); // i.e. update the [i][j] site } else { // occupied site i.e. trapping in precursor state Precursor[i][j]=1; CalculateOmega(i,j); } CalculateOmega(i,j-1); // update the nearest neighbor sites of [i][j] CalculateOmega(i,j+1); CalculateOmega(i+1,j);CalculateOmega(i-1,j); } // if K <0 end } // for j } // for i } // if K < (Depos... end } void ChooseProcess() // desorption or diffusion along the precursor state { int Flag; // = 0 found a suitable process, = 1 no process could be found int i,j; // adsorption site Flag = 1; for (i=2; ((i<=nx-2) & (Flag == 1)); i++) { // go through the lattice until a process could be found for (j=2; ((j<=ny-2) & (Flag == 1)); j++) { // go through the lattice until a process could be found if (Flag==1) K = K - Desorption[i][j]; // linear selection rule if ((K<0) & (Flag==1)) { // desorption from precursor state has been chosen as the process to proceed Flag = 0; // found a process Precursor[i][j]=0; CalculateOmega(i,j-1); // update the nearest neighbor sites of [i][j] CalculateOmega(i,j+1); CalculateOmega(i+1,j); CalculateOmega(i-1,j); } // desorption if end else { // no desorption so far if (Flag==1) K = K - Diffusion[i][j]; // linear selection rule if ((K<0) & (Flag==1)) { // diffusion event has been chosen as the process to proceed K = ran1(&idum)*(4-Neighbors[i][j]); // diffusion along precursor ok, but in what direction? // calculate new K, we have 4 possibilities // i.e. 4 neighbor sites could be free Precursor[i][j]=0; // clear current precursor site if ( (Precursor[i][j-1]==0) & (Flag==1) ) { // left hand site neighbor site free K--; // linear selection rule for diffusion process if (K<=0) { // ok, that direction has bee chosen if (Lattice[i][j-1]==0) { // diffusion hop on free lattice site Sum++; Lattice[i][j-1]=1; } else Precursor[i][j-1]=1; // adsorption site occupied so trapping in precursor CalculateOmega(i,j); // update site [i][j] CalculateOmega(i,j-1); // update the nearest neighbor sites of [i][j] CalculateOmega(i,j+1); CalculateOmega(i+1,j); CalculateOmega(i-1,j); CalculateOmega(i,j-2); // update the nearest neighbor sites of [i][j-1] CalculateOmega(i,j); CalculateOmega(i,j+1); CalculateOmega(i-1,j-1); Flag = 0; // found a process } } // if end of left hand site neighbor site free if ( (Precursor[i][j+1]==0) & (Flag==1) ) { // right hand site neighbor site free K--; if (K<=0) { if (Lattice[i][j+1]==0) { Sum++; Lattice[i][j+1]=1;} else Precursor[i][j+1]=1; CalculateOmega(i,j); CalculateOmega(i,j-1); CalculateOmega(i,j+1); CalculateOmega(i+1,j); CalculateOmega(i-1,j); CalculateOmega(i,j); CalculateOmega(i,j+2); CalculateOmega(i+1,j+1); CalculateOmega(i-1,j+1); Flag = 0; } } // right hand site if ( (Precursor[i+1][j]==0) & (Flag==1) ) { // bottom site neighbor site free K--; if (K<=0) { if (Lattice[i+1][j]==0) { Sum++; Lattice[i+1][j]=1; } else Precursor[i+1][j]=1; CalculateOmega(i,j); CalculateOmega(i,j-1); CalculateOmega(i,j+1); CalculateOmega(i+1,j); CalculateOmega(i-1,j); CalculateOmega(i,j); CalculateOmega(i-1,j); CalculateOmega(i+1,j+1); CalculateOmega(i,j+1); Flag = 0; } } // bottom hand site if ( (Precursor[i-1][j]==0) & (Flag==1) ) { // top site neighbor site free K--; if (K<=0) { if (Lattice[i-1][j]==0) { Sum++; Lattice[i-1][j]=1; } else Precursor[i-1][j]=1; CalculateOmega(i,j); CalculateOmega(i,j-1); CalculateOmega(i,j+1); CalculateOmega(i+1,j); CalculateOmega(i-1,j); CalculateOmega(i,j); CalculateOmega(i-1,j-1); CalculateOmega(i-1,j+1); CalculateOmega(i-2,j); Flag = 0; } } // top site } // diffusion if end } // else no desorption end } // for j end } // for i end } void KisliukDynamicsKMCS() { // main program segment int i,j; // loop index, adsorption site int k; // averaging loop index float CurrentCoverage; float OldCoverage; char buffer[150]; //********** initialisation Sclean = 1.0 - Pclean; Soccu = 1.0 - Poccu; ny=nx; for (k=1; k<=Average; k++) { index = 1; CurrentTime=0; OldTime=0; Sum=0; CurrentCoverage=0; OldCoverage=0; Initialize(); while (CurrentCoverage < SaturationCoverage) { Omega = Omega + Flux; dt = -log(1.0-ran1(&idum))/Omega; CurrentTime = CurrentTime + dt; K = Omega*ran1(&idum); // ran1 yields random number within ]0,1[ if (K OldCoverage)) { OldCoverage=CurrentCoverage; // time has to be elapsed even if no adsorption took place // but calculating S makes only sense if the coverage is increasing Coverage[index]=CurrentCoverage; AdsorptionProbability[index]= (nx*ny/100/Flux)/(CurrentTime - OldTime); SAverage[index]=SAverage[index]+AdsorptionProbability[index]; RealTime[index]=CurrentTime; printf("loop=%d Theta=%f t=%f S=%f SAver=%f\n", k,Coverage[index], RealTime[index], AdsorptionProbability[index],SAverage[index]); index++; OldTime=CurrentTime; }; } // while loop end } // average loop end for (i=1; i //**********How good are rundom number generators? bgh //********** parameter of the generator float U = 1; float V = 3; // 3,11,13,19,21,27,29,37 float A; // =200*U+V float B = 0; float M = 100000; float X0 = 1; // seed of the generator //********* program parameter float X; // random number int D=3; // plot points with a periode of D float XPlot,YPlot; // point coordinates int PointNumber=10; // number of points plotted float MOD(float X1, float X2) { return ((X1/X2)- (int)(X1/X2))*X2; } // modulo float Rnd() { return MOD(A*X+B,M); } // random number generator void TwoDTest() { int i,j; // loop index char answer[1]; A = 200*U+V; X = X0; for (i=1; i < PointNumber ; i++) { X = Rnd(); XPlot=X/M; for (j=1; j<=D; j++) { X = Rnd();}; YPlot=X/M; printf("%f %f \n",XPlot,YPlot);} scanf("%s",&answer); } void main () { TwoDTest(); } **************************************************************************************************** //////////////////////////////////////////////////////////////////////////////////////////////////// **************************************************************************************************** Listing A7: #include // bgh #include #include #define NRANSI #include "write.h" int nx = 9; // size of the lattice int ny = 9; // nx=ny=const=9 only for the demo int List[100]; // list structure int TwoDLattice[10][10]; // corresponding matrix int XR1[101]; // neigbor list x-directing, right int XL1[101]; // x-dir. left int YL1[101]; // y-dir. left int YR1[101]; // y-dir. right char answer[1]; // for simple screen output int i,j,k; // loop index void PeriodicBoundary() { long int I; for (I = 1; I <= nx*ny; I++) { XR1[I] = I + 1; XL1[I] = I - 1; YL1[I] = I - nx; YR1[I] = I + nx; } for (I = 1; I <= nx; I++) { YL1[I] = nx * (ny - 1) + I; YR1[nx * (ny - 1) + I] = I; } for (I = 1; I <= ny; I++) { XL1[1 + nx * (I - 1)] = nx * I; XR1[nx * I] = nx * (I - 1) + 1; } } void InitData () { // set all lists and arrays to zero for (i=1; i<=(nx*ny); i++) {List[i]= 0;} for (i=1; i<=nx; i++) { for (j=1; j<=ny; j++) { TwoDLattice[i][j]=0; }; }; } void GenerateList () { // generate list such that if fullfills the nomenclature typically used for a matrix for (i=1; i<=nx; i++) { for (j=1; j<=ny; j++) { List[k]= i*10+j; k++; } } } void ListToMatrix () { // convert list to matrix, typically required for plotting a configuration k=1; for (i=1; i<=nx; i++) { for (j=1; j<=ny; j++) { TwoDLattice[i][j]=List[nx*(i-1)+j]; } } } void ShowListAsAList () { for (i=1; i<=(nx*ny); i++) { printf("%d ",List[i]);} } void ShowListAsAMatrix () { // this screen output works for (i=1; i<=(nx*ny); i++) { // only for a 9 x 9 matrix printf("%d ",List[i]); if ( nx*i % (nx*ny) == 0 ) printf("\n"); } } void ShowMatrix () { for (i=1; i<=nx; i++) { for (j=1; j<=ny; j++) { printf("[%d,%d]=%d ",i,j,TwoDLattice[i][j]); } printf("\n"); } } void ListDemo() { int I; // number of neigbor list element printf("List demo\n\n"); InitData(); GenerateList(); printf("Show the list \n"); ShowListAsAList(); printf("\nPress a key & return\n");scanf("%s",answer); printf("\nShow the list as a matrix \n"); printf("11: row = 1 column = 1\n"); printf("32: row = 3 column = 2\n\n"); ShowListAsAMatrix(); printf("\n"); printf("\nPress a key & return\n");scanf("%s",answer); printf("Convert list to matrix and show matrix \n"); ListToMatrix(); ShowMatrix(); printf("\nPress a key & return\n");scanf("%s",answer); PeriodicBoundary(); I = 1; while (I!=999) { printf("Demo of neigbor lists \n (finish demo? Then type 999)\n"); printf("Chose a list element number smaller than %d = ",nx*ny); scanf("%d",&I); if (I!=999) { printf("\nList no. = %d element [%d]=%d\n",I,I,List[I]); printf("x-right = %d x-left = %d y-left = %d y-right = %d\n", List[XR1[I]], List[XL1[I]], List[YL1[I]], List[YR1[I]]); printf("Correct?\n"); ShowListAsAMatrix(); printf("It should be correct\n"); } else printf("By by \n"); } } void main() { ListDemo(); } **************************************************************************************************** //////////////////////////////////////////////////////////////////////////////////////////////////// **************************************************************************************************** Listing 8: int random = 2 * 5731 - 1; // Ste float RND() { // random number equally distributed between 0 and 1 random = random * 16807; Var = random * 2.328306E-10 + 0.5; return Var; } float GAUSS() { // Gauss distributed random numbers float y1, y2, s; s = 1.0; while (s >=1.0) { y1 = (RND()-0.5)*2; y2 = (RND()-0.5)*2; s = (y1 * y1) + (y2 * y2); } return y1 * (sqrt) (- (log)(s) / s) ; } double RAND0 = 0.12345678901; // Rogo double Random() // random number within [0,1) { RAND0 = 3125.0 * RAND0; RAND0 = RAND0 - (double)((int)RAND0); return ( RAND0 ); } **************************************************************************************************** //////////////////////////////////////////////////////////////////////////////////////////////////// **************************************************************************************************** Listing : // simple demo of the random number generator from the // numerical recipes library bgh //e.g. the ran1 function generates random numbers within the interval (0,1) int idum=(-1); // for the random number generator, // include always this line in the examples void RanDemo() { int i; int X,Y; for (i=1; (i < 100) ; i++) { X = ran1(&idum)*100; Y = ran4(&idum)*100; printf("%d %d \n",X,Y); } } main (void) { RanDemo();} **************************************************************************************************** //////////////////////////////////////////////////////////////////////////////////////////////////// **************************************************************************************************** Listing : #include #include #include // standard "C" header #include // for simple mathematical functions #include // for simple input and output operations #include "nr.h" // from numerical recipes #define NRANSI // some compiler require this statement, // which indicates that standard "C" code follows // simple demo of the random number generator from the // numerical recipes library //e.g. the ran1 function generates random numbers within the interval (0,1) int idum=(-1); // for the random number generator, // include always this line in the examples int random = 2 * 5731 - 1; // Ste float RND() { // random number equally distributed between 0 and 1 float Var; random = random * 16807; Var = random * 2.328306E-10 + 0.5; return Var; } float GAUSS() { // Gauss distributed random numbers float y1, y2, s; s = 1.0; while (s >=1.0) { y1 = (RND()-0.5)*2; y2 = (RND()-0.5)*2; s = (y1 * y1) + (y2 * y2); } return y1 * (sqrt) (- (log)(s) / s) ; } double RAND0 = 0.12345678901; // Rogo double Random() // random number within [0,1) { RAND0 = 3125.0 * RAND0; RAND0 = RAND0 - (double)((int)RAND0); return ( RAND0 ); } void RanDemo() { int i; int X,Y; for (i=1; (i < 10) ; i++) { X = ran4(&idum)*100; Y = ran4(&idum)*100; printf("%d %d %f %f\n",X,Y,Random(),RND()); } } main (void) { RanDemo();} **************************************************************************************************** //////////////////////////////////////////////////////////////////////////////////////////////////// **************************************************************************************************** Dr. Uwe Burghaus Assistant Professor - Surface Chemistry North Dakota State University, Department of Chemistry and Molecular Biology Fargo, North Dakota 58105-5516 Phone 701-231-9742 FAX 701-231-8831 Email Uwe.Burghaus@ndsu.edu http://www.chem.ndsu.nodak.edu http://www.uweburghaus.de