/* Steven Andrews, started 10/22/01, many modifications since then. This is the library of routines for the Smoldyn program */ /* See documentation called Smoldyn_doc1.doc and Smoldyn_doc2.doc */ /* Copyright 2003-2006 by Steven Andrews. Permission is granted for non-commercial use of and modifications to the code. */ #include #include #include #include #include #include "Rn.h" #include "Zn.h" #include "random.h" #include "string2.h" #include "math2.h" #include "rxnparam.h" #include "SimCommand.h" #include "smollib.h" #include "smollib2.h" #include "VoidComp.h" #include "Geometry.h" #include "opengl2.h" // See the documentation for a discussion of macros and global variables. #define RANDTABLEMAX 4095 #define CHECK(A) if(!(A)) goto failure #define CHECKS(A,B) if(!(A)) {strncpy(erstr,B,STRCHAR);goto failure;} double GaussTable[RANDTABLEMAX+1]; /******************************************************************************/ /********************************* Molecules **********************************/ /******************************************************************************/ /* molalloc allocates and initiallizes a new moleculestruct. The box and diffusion matrix members are returned as NULL. */ moleculeptr molalloc(int dim) { moleculeptr mptr; int d; if(dim<1) return NULL; mptr=(moleculeptr) malloc(sizeof(struct moleculestruct)); if(!mptr) return NULL; mptr->serno=0; mptr->pos=NULL; mptr->posx=NULL; mptr->wrap=NULL; mptr->ident=0; mptr->box=NULL; mptr->pnl=NULL; mptr->face='f'; mptr->via=NULL; mptr->pos=(double*)calloc(dim,sizeof(double)); if(!mptr->pos) {molfree(mptr);return NULL;} mptr->posx=(double*)calloc(dim,sizeof(double)); if(!mptr->posx) {molfree(mptr);return NULL;} mptr->wrap=(int*)calloc(dim,sizeof(int)); if(!mptr->wrap) {molfree(mptr);return NULL;} mptr->via=(double*)calloc(dim,sizeof(double)); if(!mptr->via) {molfree(mptr);return NULL;} for(d=0;dpos[d]=mptr->posx[d]=mptr->via[d]=0; mptr->wrap[d]=0; } return mptr; } /* molfree frees the space allocated for a moleculestruct, as well as its two position vectors. */ void molfree(moleculeptr mptr) { if(!mptr) return; if(mptr->pos) free(mptr->pos); if(mptr->posx) free(mptr->posx); if(mptr->wrap) free(mptr->wrap); if(mptr->via) free(mptr->via); free(mptr); return; } /* molssalloc allocates and initiallizes a molecule superstructure with max molecule spaces in each of the three lists. max must be at least 1. The dead list is filled with empty molecules. The molecule boxes are left as NULLs, and need to be set. Book keeping elements for the lists are set to their initial values. */ molssptr molssalloc(int dim,int max,int maxident) { molssptr mols; int i,m,ll; if(max<1) return NULL; mols=(molssptr) malloc(sizeof(struct molsuperstruct)); if(!mols) return NULL; mols->difc=NULL; mols->difstep=NULL; mols->difm=NULL; mols->display=NULL; mols->color=NULL; mols->grphdata=NULL; mols->grphser=NULL; mols->live[0]=NULL; mols->live[1]=NULL; mols->dead=NULL; mols->max=max; mols->nl[0]=0; mols->nl[1]=0; mols->topl[0]=0; mols->topl[1]=0; mols->nd=max; mols->topd=max; mols->serno=1; CHECK(mols->difc=(double*) calloc(maxident,sizeof(double))); CHECK(mols->difstep=(double*) calloc(maxident,sizeof(double))); for(i=0;idifc[i]=mols->difstep[i]=0; CHECK(mols->difm=(double **) calloc(maxident,sizeof(double *))); for(i=0;idifm[i]=NULL; CHECK(mols->display=(double*) calloc(maxident,sizeof(double))); for(i=0;idisplay[i]=3.0; CHECK(mols->color=(double **) calloc(maxident,sizeof(double *))); for(i=0;icolor[i]=NULL; for(i=0;icolor[i]=(double*) calloc(3,sizeof(double))); mols->color[i][0]=mols->color[i][1]=mols->color[i][2]=0; } CHECK(mols->grphdata=(double **) calloc(max,sizeof(double *))); for(m=0;mgrphdata[m]=NULL; for(m=0;mgrphdata[m]=(double*) calloc(3,sizeof(double))); mols->grphdata[m][0]=mols->grphdata[m][1]=mols->grphdata[m][2]=0; } CHECK(mols->grphser=(int*) calloc(max,sizeof(int))); for(m=0;mgrphser[m]=0; for(ll=0;ll<2;ll++) { CHECK(mols->live[ll]=(moleculeptr *) calloc(max,sizeof(moleculeptr))); for(m=0;mlive[ll][m]=NULL; } CHECK(mols->dead=(moleculeptr *) calloc(max,sizeof(moleculeptr))); for(m=0;mdead[m]=NULL; for(m=0;mdead[m]=molalloc(dim)); return mols; failure: molssfree(mols,maxident); return NULL; } /* molssfree frees both a superstructure of molecules and all the molecules in all its lists. */ void molssfree(molssptr mols,int maxident) { int m,ll,i; if(!mols) return; if(mols->difc) free(mols->difc); if(mols->difstep) free(mols->difstep); if(mols->difm) { for(i=0;idifm[i]) freeM(mols->difm[i]); free(mols->difm); } if(mols->display) free(mols->display); if(mols->color) { for(i=0;icolor[i]) free(mols->color[i]); free(mols->color); } if(mols->grphdata) { for(m=0;mmax;m++) if(mols->grphdata[m]) free(mols->grphdata[m]); free(mols->grphdata); } if(mols->grphser) free(mols->grphser); for(ll=0;ll<2;ll++) if(mols->live[ll]) { for(m=0;mnl[ll];m++) if(mols->live[ll][m]) molfree(mols->live[ll][m]); free(mols->live[ll]); } if(mols->dead) { for(m=0;mnd;m++) if(mols->dead[m]) molfree(mols->dead[m]); free(mols->dead); } free(mols); return; } /* molssoutput prints all the parameters in a molecule superstructure. While it should not be needed, hopefully, this routine looks for and prints out information on molecules that are not sorted correctly in the live and dead lists. It also prints out information about each molecule, including diffusion constants, rms step lengths, colors, and display sizes. */ void molssoutput(simptr sim) { int nident,d,m,ll,i,outside; molssptr mols; char **name; printf("MOLECULE PARAMETERS\n"); if(!sim||!sim->mols) { printf(" No molecule superstructure defined\n\n"); return; } mols=sim->mols; name=sim->name; nident=sim->nident; printf(" %i molecules allocated\n",mols->max); printf(" Live mobile: %i, live fixed: %i\n",mols->nl[0],mols->nl[1]); printf(" Dead: %i",mols->nd); if(mols->topd==mols->nd) printf("\n"); else printf(" Top of empty molecule list: %i\n",mols->topd); for(ll=0;ll<2;ll++) { for(m=0;mnl[ll];m++) if(!mols->live[ll][m]) printf(" BUG: NULL molecule in live list %i at %i\n",ll,m); else if(!mols->live[ll][m]->ident) printf(" WARNING: empty molecule in live list %i at %i\n",ll,m); for(;mmax;m++) if(mols->live[ll][m]) printf(" BUG: misplaced molecule in live list %i at %i\n",ll,m); } for(m=0;mtopd;m++) if(!mols->dead[m]) printf(" BUG: NULL molecule in dead list at %i\n",m); else if(mols->dead[m]->ident) printf(" BUG: live molecule in dead list at %i\n",m); for(;mnd;m++) if(!mols->dead[m]) printf(" BUG: NULL molecule in resurrected list at %i\n",m); else if(!mols->dead[m]->ident) printf(" BUG: dead molecule in resurrected list at %i\n",m); for(;mmax;m++) if(mols->dead[m]) printf(" BUG: misplaced molecule in dead list at %i\n",m); outside=0; for(ll=0;ll<2;ll++) for(m=0;mnl[ll];m++) for(d=0;ddim;d++) if(mols->live[ll][m]->pos[d]wlist[2*d]->pos||mols->live[ll][m]->pos[d]>sim->wlist[2*d+1]->pos) outside++; if(outside) printf(" WARNING: %i molecules are outside system walls\n",outside); for(i=1;idifc[i],mols->difstep[i],mols->difm[i]?"anisotropic":" "); if(sim->graphics&&mols->display[i]) printf(" color= %lf %lf %lf, size= %lf\n",mols->color[i][0],mols->color[i][1],mols->color[i][2],mols->display[i]); else if(sim->graphics) printf(" not displayed to graphics\n"); } printf("\n"); return; } /* This sets up the diffusion coefficient related aspects of the molecule superstructure. If difm is initiallized but not difc, then it creates difc. It also assigns difstep. */ void setdiffusion(simptr sim) { int i,nident,dim; molssptr mols; nident=sim->nident; dim=sim->dim; mols=sim->mols; for(i=0;idifm[i]&&!mols->difc[i]) { dotMMD(mols->difm[i],mols->difm[i],sim->m2,dim,dim,dim); mols->difc[i]=traceMD(sim->m2,dim)/dim; } for(i=0;idifstep[i]=sqrt(2.0*mols->difc[i]*sim->dt); return; } /* molsort updates the live and dead lists of both a molecule superstructure and the relevent boxes after a reaction or other changes. If difsort is 0, it assumes that live molecules are in the correct live list, otherwise, it sorts this too (the only way itÕs likely to need sorting is if a command changes a moleculeÕs identity). First it might deal with diffusion sorting by moving missorted molecules to the resurrected list. Afterwards, and in normal operation, it moves the resurrected molecules off the top of the dead list (those numbered between topd and nd-1) to the live list. Then it moves the expired molecules from the live list to the dead list. Finally it compacts the live list. Molecule ordering in lists is not preserved. It is required that the box numbers for molecules are correct and that the molecule is listed in the box list corresponding to the master live lists. The routine returns 0 for normal operation and 1 if memory could not be allocated when a box was being expanded. */ int molsort(molssptr mols,int difsort) { int m,m2,m3,ll,reborn[2]; moleculeptr mptr,*dead,*live; boxptr bptr; double *difc; int *nl; difc=mols->difc; dead=mols->dead; nl=mols->nl; if(difsort) // sort molecules missorted by diffusion onto resurrected list for(ll=0;ll<2;ll++) { live=mols->live[ll]; for(m=0;mident]==0)) { dead[mols->nd++]=live[m]; live[m]=NULL; }} reborn[0]=reborn[1]=0; for(m=mols->topd;mnd;m++) { // move resurrected molecules onto live lists mptr=dead[m]; ll=mols->difc[mptr->ident]==0; mols->live[ll][mols->nl[ll]++]=mptr; mols->dead[m]=NULL; mptr->serno=mols->serno++; reborn[ll]++; bptr=mptr->box; if(bptr->nmol[ll]==bptr->maxmol[ll]) if(expandbox(bptr,5,ll)) return 1; bptr->mol[ll][bptr->nmol[ll]++]=mptr; } mols->nd=mols->topd; for(ll=0;ll<2;ll++) { // move expired molecules from live to dead list live=mols->live[ll]; m2=nl[ll]; for(m=nl[ll]-1;m>=0;m--) if(!live[m]->ident) { mptr=live[m]; dead[mols->nd++]=mptr; live[m2=m]=NULL; bptr=mptr->box; for(m3=bptr->nmol[ll]-1;bptr->mol[ll][m3]!=mptr;m3--); bptr->mol[ll][m3]=bptr->mol[ll][--bptr->nmol[ll]]; } for(m=m2++;mtopl[ll]=nl[ll]-reborn[ll]; } mols->topd=mols->nd; return 0; } /* diffuse does the diffusion for all molecules in the live[0] list (mobile), over one time step. Walls are ignored and molecules are not reassigned to the boxes. If there is a diffusion matrix, it is used for anisotropic diffusion; otherwise isotropic diffusion is done, using the difstep parameter. */ void diffuse(simptr sim) { int m,d,nmol,dim,i; double flt1; double *v1,*v2,*difstep,**difm; moleculeptr *mlist; dim=sim->dim; v1=sim->v1; v2=sim->v2; mlist=sim->mols->live[0]; nmol=sim->mols->nl[0]; difstep=sim->mols->difstep; difm=sim->mols->difm; flt1=sqrt(2.0*sim->dt); for(m=0;mident; if(!difm[i]) for(d=0;dposx[d]=mlist[m]->pos[d]; mlist[m]->pos[d]+=difstep[i]*GaussTable[rand()&RANDTABLEMAX]; } else { for(d=0;dposx[d]=mlist[m]->pos[d]; v1[d]=flt1*GaussTable[rand()&RANDTABLEMAX]; } dotMVD(difm[i],v1,v2,dim,dim); for(d=0;dpos[d]+=v2[d]; }} return; } /******************************************************************************/ /************************************ Walls ***********************************/ /******************************************************************************/ /* wallalloc allocates and initializes a new wall. The pointer to the opposite wall needs to be set. */ wallptr wallalloc(void) { wallptr wptr; wptr=(wallptr) malloc(sizeof(struct wallstruct)); if(!wptr) return 0; wptr->wdim=0; wptr->side=0; wptr->pos=0; wptr->type='r'; wptr->opp=NULL; return wptr; } /* wallfree frees a wall. */ void wallfree(wallptr wptr) { if(!wptr) return; free(wptr); return; } /* wallsalloc allocates an array of pointers to 2*dim walls, allocates each of the walls, and sets them to default conditions (reflecting walls at 0 and 1 on each coordinate) with correct pointers in each opp member. */ wallptr *wallsalloc(int dim) { int w,d; wallptr *wlist; if(dim<1) return NULL; wlist=(wallptr *) calloc(2*dim,sizeof(wallptr)); if(!wlist) return 0; for(w=0;w<2*dim;w++) wlist[w]=NULL; for(w=0;w<2*dim;w++) if(!(wlist[w]=wallalloc())) {wallsfree(wlist,dim);return 0;} for(d=0;dwdim=wlist[2*d+1]->wdim=d; wlist[2*d]->side=0; wlist[2*d+1]->side=1; wlist[2*d]->pos=0; wlist[2*d+1]->pos=1; wlist[2*d]->type=wlist[2*d+1]->type='r'; wlist[2*d]->opp=wlist[2*d+1]; wlist[2*d+1]->opp=wlist[2*d]; } return wlist; } /* wallsfree frees an array of 2*dim walls, including the walls. */ void wallsfree(wallptr *wlist,int dim) { if(!wlist||dim<1) return; for(dim--;dim>=0;dim--) { wallfree(wlist[2*dim+1]); wallfree(wlist[2*dim]); } return; } /* walloutput prints the wall structure information, including wall dimensions, positions, and types, as well as the total simulation volume. */ void walloutput(int dim,wallptr *wlist) { int w,d; wallptr wptr; double vol; printf("WALL PARAMETERS\n"); if(!wlist) { printf(" No walls defined for simulation\n\n"); return; } for(w=0;w<2*dim;w++) { wptr=wlist[w]; printf(" Wall %i: dimension %i, at %lf, type %c\n",w,wptr->wdim,wptr->pos,wptr->type); if(wlist[w+1-2*(w%2)]!=wptr->opp) printf(" ERROR: opposing wall is incorrect\n"); } vol=1; for(d=0;dpos-wlist[2*d]->pos; if(dim==1) printf(" Length: %lf\n",vol); else if(dim==2) printf(" Area: %lf\n",vol); else printf(" Volume: %lf\n",vol); printf("\n"); return; } /* checkwalls does the reflection, wrap-around, or absorption of molecules at walls by checking the current position, relative to the wall positions (as well as a past position for absorbing walls). It does not reassign the molecules to boxes or sort the live and dead ones, although this typically will not be required anyhow. It returns the number of wall collisions that were detected and processed during that time step. */ int checkwalls(simptr sim) { int w,b,wctr,nmol,nbox,m,d; boxptr *blist; moleculeptr *mlist; wallptr wptr; double pos2,diff,difi,step,*difstep; if(sim->srfss) return 0; wctr=0; nbox=sim->boxs->nbox; blist=sim->boxs->blist; difstep=sim->mols->difstep; for(b=0;bnwall&&blist[b]->nmol[0]) { nmol=blist[b]->nmol[0]; mlist=blist[b]->mol[0]; for(w=0;wnwall;w++) { wptr=blist[b]->wlist[w]; d=wptr->wdim; if(wptr->type=='r'&&!(wptr->side)) { // reflective pos2=2*wptr->pos; for(m=0;mpos[d]pos) { wctr++; mlist[m]->pos[d]=pos2-mlist[m]->pos[d];}} else if(wptr->type=='r') { pos2=2*wptr->pos; for(m=0;mpos[d]>wptr->pos) { wctr++; mlist[m]->pos[d]=pos2-mlist[m]->pos[d];}} else if(wptr->type=='p'&&!(wptr->side)) { // periodic pos2=wptr->opp->pos-wptr->pos; for(m=0;mpos[d]pos) { wctr++; mlist[m]->pos[d]=pos2+mlist[m]->pos[d]; mlist[m]->wrap[d]--; }} else if(wptr->type=='p') { pos2=wptr->opp->pos-wptr->pos; for(m=0;mpos[d]>wptr->pos) { wctr++; mlist[m]->pos[d]=pos2+mlist[m]->pos[d]; mlist[m]->wrap[d]++; }} else if(wptr->type=='a') { // absorbing for(m=0;mpos-mlist[m]->pos[d]; difi=wptr->pos-mlist[m]->posx[d]; step=difstep[mlist[m]->ident]; if((!(wptr->side)&&diff>0)||(wptr->side&&diff<0)||unirand30(0,1)ident=0; }}}}} return wctr; } /******************************************************************************/ /********************************** Reactions *********************************/ /******************************************************************************/ /* rxnalloc allocates and initializes a reaction structure, leaving it fully set up but with zero reactions. It can be used as is, but it wonÕt do anything until reactions are added. order and total need to be at least 0 and maxident needs to be at least 1. */ rxnptr rxnalloc(int order,int maxident,int total) { rxnptr rxn; int i,r,ni2o; rxn=NULL; CHECK(order>=0&&maxident>0&&total>=0); CHECK(rxn=(rxnptr) malloc(sizeof(struct rxnstruct))); rxn->order=order; rxn->nrxn=NULL; rxn->table=NULL; rxn->lists=0; rxn->total=total; rxn->rname=NULL; rxn->rate=NULL; rxn->rate2=NULL; rxn->rpar=NULL; rxn->rpart=NULL; rxn->nprod=NULL; rxn->prod=NULL; ni2o=intpower(maxident,order); CHECK(rxn->nrxn=(int*)calloc(ni2o,sizeof(int))); for(i=0;inrxn[i]=0; CHECK(rxn->table=(int **) calloc(ni2o,sizeof(int *))); for(i=0;itable[i]=NULL; if(total) { CHECK(rxn->rname=(char **)calloc(total,sizeof(char *))); for(r=0;rrname[r]=NULL; for(r=0;rrname[r]=EmptyString()); CHECK(rxn->rate=(double*) calloc(total,sizeof(double))); CHECK(rxn->rate2=(double*) calloc(total,sizeof(double))); CHECK(rxn->rpar=(double*) calloc(total,sizeof(double))); CHECK(rxn->rpart=(char *) calloc(total,sizeof(char))); for(r=0;rrate[r]=rxn->rate2[r]=-1; for(r=0;rrpar[r]=0; for(r=0;rrpart[r]='?'; CHECK(rxn->nprod=(int*)calloc(total,sizeof(int))); for(r=0;rnprod[r]=0; CHECK(rxn->prod=(moleculeptr **) calloc(total,sizeof(moleculeptr *))); for(r=0;rprod[r]=NULL; } return rxn; failure: rxnfree(rxn,maxident); return NULL; } /* rxnfree frees a reaction structure for any order reaction. */ void rxnfree(rxnptr rxn,int maxident) { int r,p,i; int order,total,*nrxn,**table,*nprod; moleculeptr **prod; if(!rxn) return; order=rxn->order; total=rxn->total; nrxn=rxn->nrxn; table=rxn->table; nprod=rxn->nprod; prod=rxn->prod; if(nprod&&prod) for(r=0;rrname) { for(r=0;rrname[r]) free(rxn->rname[r]); free(rxn->rname); } if(rxn->rate) free(rxn->rate); if(rxn->rate2) free(rxn->rate2); if(rxn->rpar) free(rxn->rpar); if(rxn->rpart) free(rxn->rpart); if(nrxn&&table) for(i=0;idim; maxident=sim->maxident; nident=sim->nident; name=sim->name; while(fgets(line,STRCHAR,fptr)) { (*lctrptr)++; if(strchr(line,'#')) *strchr(line,'#')='\0'; itct=sscanf(line,"%s",word); if(itct) line2=strnword(line,2); else line2=NULL; if(itct<=0); // blank line else if(!strcmp(word,"/*")) { // start comment with /* incomment=1; } else if(!strcmp(word,"*/")) { // end comment with */ CHECKS(incomment,"*/ without a corresponding /*"); incomment=0; } else if(incomment); else if(!strcmp(word,"end_reaction")) { // end_reaction CHECKS(!line2,"unexpected text following end_reaction"); return 0; } else if(!line2) { // just word CHECKS(0,"unknown word or missing parameter"); } else if(!strcmp(word,"order")) { // order, got[0] CHECKS(!got[0],"order can only be entered once"); got[0]++; itct=sscanf(line2,"%i",&order); CHECKS(itct==1,"error reading order"); CHECKS(order>=0&&order<=2,"order needs to be between 0 and 2"); if(sim->rxn[order]) { got[1]++; rxn=sim->rxn[order]; nrxn=rxn->nrxn; table=rxn->table; rname=rxn->rname; rate=rxn->rate; rate2=rxn->rate2; rpar=rxn->rpar; rpart=rxn->rpart; nprod=rxn->nprod; prod=rxn->prod; } CHECKS(!strnword(line2,2),"unexpected text following order"); } else if(!strcmp(word,"max_rxn")) { // max_rxn, got[1] CHECKS(got[0],"order needs to be entered before max_rxn"); CHECKS(!got[1],"max_rxn can only be entered once"); got[1]++; itct=sscanf(line2,"%i",&total); CHECKS(itct==1,"error reading max_rxn"); CHECKS(total>=0,"num_rxn needs to be >=0"); sim->rxn[order]=rxn=rxnalloc(order,maxident,total); CHECKS(rxn,"memory allocation error in rxnalloc()"); nrxn=rxn->nrxn; table=rxn->table; rname=rxn->rname; rate=rxn->rate; rate2=rxn->rate2; rpar=rxn->rpar; rpart=rxn->rpart; nprod=rxn->nprod; prod=rxn->prod; CHECKS(!strnword(line2,2),"unexpected text following max_rxn"); } else if(!strcmp(word,"reactant")&&order==0) { // reactant, 0 CHECKS(got[1],"max_rxn needs to be entered before reactant"); j=wordcount(line2); CHECKS(j>0,"number of reactions needs to be >0"); for(j--;j>=0;j--) { itct=sscanf(line2,"%s",nm); CHECKS(itct==1,"missing reaction name in reactant"); r=stringfind(rname,total,nm); if(r<0) { r=stringfind(rname,total,""); CHECKS(r>=0,"not enough zeroth order reactions allocated"); strncpy(rname[r],nm,STRCHAR); } if(j) CHECKS(line2=strnword(line2,2),"reaction names missing in reactant"); }} else if(!strcmp(word,"reactant")&&order==1) { // reactant, 1 CHECKS(got[1],"max_rxn needs to be entered before reactant"); itct=sscanf(line2,"%s",nm); CHECKS(itct==1,"format for reactant: mol_name rxn_list"); i=stringfind(name,nident,nm); CHECKS(i>0,"name of molecule not recognized for reactant"); j=wordcount(line2)-1; CHECKS(j>0,"number of reactions needs to be >0"); iptr=table[i]; j1=nrxn[i]; nrxn[i]+=j; CHECKS(table[i]=allocZV(nrxn[i]),"memory allocation failed in allocZV()"); for(j=0;j=0,"not enough first order reactions allocated"); strncpy(rname[r],nm,STRCHAR); } table[i][j]=r; }} else if(!strcmp(word,"reactant")&&order==2) { // reactant, 2 CHECKS(got[1],"max_rxn needs to be entered before reactants"); itct=sscanf(line2,"%s + %s",nm,nm2); CHECKS(itct==2,"format for reactant: name + name rxn_list"); i1=stringfind(name,nident,nm); i2=stringfind(name,nident,nm2); CHECKS(i1>0&&i2>0,"name not recognized for reactant"); j=wordcount(line2)-3; CHECKS(j>0,"number of reactions needs to be >0"); i=i1*maxident+i2; iptr=table[i]; j1=nrxn[i]; nrxn[i]+=j; CHECKS(table[i]=allocZV(nrxn[i]),"memory allocation failed in allocZV()"); for(j=0;j=0,"not enough second order reactions allocated"); strncpy(rname[r],nm,STRCHAR); } table[i][j]=r; } if(i1!=i2) { i=i2*maxident+i1; nrxn[i]=nrxn[i1*maxident+i2]; if(table[i]) freeZV(table[i]); CHECKS(table[i]=allocZV(nrxn[i]),"memory allocation failed in allocZV()"); copyZV(table[i1*maxident+i2],table[i],nrxn[i]); }} else if(!strcmp(word,"rate")) { // rate CHECKS(got[1],"max_rxn needs to be entered before rate"); itct=sscanf(line2,"%s %lf",nm,&rtemp); CHECKS(itct==2,"format for rate: rxn_name rate"); r=stringfind(rname,total,nm); CHECKS(r>=0,"unknown reaction name in rate"); rate[r]=rtemp; CHECKS(rate[r]>=0,"reaction rate needs to be >=0 (maybe try rate_internal)"); CHECKS(!strnword(line2,3),"unexpected text following rate"); } else if(!strcmp(word,"rate_internal")) { // rate_internal CHECKS(got[1],"max_rxn needs to be entered before rate_internal"); itct=sscanf(line2,"%s %lf",nm,&rtemp); CHECKS(itct==2,"format for rate_internal: rxn_name rate"); r=stringfind(rname,total,nm); CHECKS(r>=0,"unknown reaction name in rate"); rate2[r]=rtemp; CHECKS(rate2[r]>=0,"rate_internal needs to be >=0"); CHECKS(!strnword(line2,3),"unexpected text following rate_internal"); } else if(!strcmp(word,"product")) { // product CHECKS(got[1],"max_rxn needs to be entered before product"); itct=sscanf(line2,"%s",nm); CHECKS(itct==1,"format for product: rxn_name product_list"); r=stringfind(rname,total,nm); CHECKS(r>=0,"unknown reaction name in product"); nptemp=wordcount(line2)/2; CHECKS(nptemp>=0,"number of products needs to be >=0"); CHECKS(nprod[r]==0,"products for a reaction can only be entered once"); nprod[r]=nptemp; if(nprod[r]) { prod[r]=(moleculeptr*) calloc(nprod[r],sizeof(moleculeptr)); CHECKS(prod[r],"memory allocation failed for product"); for(p=0;p0,"name not recognized for product"); mptr->ident=i; if(p+1=0,"unknown reaction name in product_param"); CHECKS(line2=strnword(line2,2),"format for product_param: rxn type [parameters]"); itct=sscanf(line2,"%c",&rpart[r]); CHECKS(itct==1&&strchr("pxrbqysofi",rpart[r]),"invalid parameter type in product_param"); if(strchr("pxrbqys",rpart[r])) { CHECKS(line2=strnword(line2,2),"missing parameter in product_param"); itct=sscanf(line2,"%lf",&rpar[r]); CHECKS(itct==1,"error reading parameter in product_param"); if(rpart[r]=='p'||rpart[r]=='q') { CHECKS(rpar[r]>0&&rpar[r]<1,"product_param value needs to be between 0 and 1"); } else if(rpart[r]=='x'||rpart[r]=='y') { CHECKS(rpar[r]>0&&rpar[r]<=1,"product_param value needs to be between 0 and 1"); } else CHECKS(rpar[r]>=0,"product_param value needs to be 0 or larger"); CHECKS(!strnword(line2,2),"unexpected text following product_param"); } else if(strchr("of",rpart[r])) { CHECKS(line2=strnword(line2,2),"missing parameters in product_param"); itct=sscanf(line2,"%s",nm2); CHECKS(itct==1,"format for product_param: rxn type [parameters]"); CHECKS((i=stringfind(name,nident,nm2))>=0,"unknown molecule in product_param"); for(p=0;pident!=i;p++); CHECKS(ppos,NULL); CHECKS(itct==dim,"insufficient data for position vector for product_param"); CHECKS(!strnword(line2,dim+1),"unexpected text following product_param"); } else { CHECKS(!strnword(line2,2),"unexpected text following product_param"); }} else { // unknown word CHECKS(0,"syntax error within reaction definition"); }} CHECKS(0,"end of file encountered before end statement"); // end of file failure: // failure rxnfree(rxn,maxident); sim->rxn[order]=NULL; return *lctrptr; } /* rxnoutput displays the complete contents of a reaction structure for order order. It also does some other calculations, such as the probability of geminate reactions for the products and the diffusion and activation limited rate constants. */ void rxnoutput(simptr sim,int order) { int i,j,j3,r,r3,p,i1,i2,i3,maxident,nident,rev; int dim,total,*nrxn,**table,*nprod; double *rate,*rate2,*rpar,*difc,dsum,step,bval,pgem; moleculeptr **prod; char *rpart,**name,**rname; rxnptr rxn; printf("ORDER %i REACTION PARAMETERS\n",order); if(!sim||!sim->rxn[order]) { printf(" No reactions of order %i\n\n",order); return; } rxn=sim->rxn[order]; dim=sim->dim; maxident=sim->maxident; nident=sim->nident; name=sim->name; order=rxn->order; total=rxn->total; nrxn=rxn->nrxn; table=rxn->table; nprod=rxn->nprod; rname=rxn->rname; rate=rxn->rate; rate2=rxn->rate2; rpar=rxn->rpar; rpart=rxn->rpart; prod=rxn->prod; difc=sim->mols->difc; printf(" %i reactions allocated\n",total); if(order==0) { for(r=0;r ",rname[r]); for(p=0;pident],p+1pos,prod[r][p]->pos,dim)>0) { printf(" product %s at ",name[prod[r][p]->ident]); printVD(prod[r][p]->pos,dim); } if(rate[r]>=0) printf(" requested rate constant: %lf\n",rate[r]); if(rpart[r]!='?') printf(" product type, parameter: %c %lf\n",rpart[r],rpar[r]); printf(" average reactions per time step: %lf\n",rate2[r]); printf(" actual rate constant: %lf\n",calcrate(sim,0,0,r,NULL)); if(nprod[r]==2&&sim->rxn[2]&&sim->rxn[2]->nrxn[i3=prod[r][0]->ident*maxident+prod[r][1]->ident]) for(j3=0;j3rxn[2]->nrxn[i3];j3++) { r3=sim->rxn[2]->table[i3][j3]; if(sim->rxn[2]->rate2[r3]>=0) { bval=distanceVVD(prod[r][0]->pos,prod[r][1]->pos,dim); dsum=sim->mols->difc[prod[r][0]->ident]+sim->mols->difc[prod[r][1]->ident]; step=sqrt(2.0*sim->dt*dsum); pgem=1.0-numrxnrate(step,sqrt(sim->rxn[2]->rate2[r3]),-1)/numrxnrate(step,sqrt(sim->rxn[2]->rate2[r3]),bval); printf(" prob. of geminate continuation reaction '%s' is %lf\n",sim->rxn[2]->rname[r3],pgem); }}}} else if(order==1) { for(i=0;i ",rname[r],name[i]); for(p=0;pident],p+1pos,prod[r][p]->pos,dim)>0) { printf(" product %s at ",name[prod[r][p]->ident]); printVD(prod[r][p]->pos,dim); } if(rate[r]>=0) printf(" requested rate constant: %lf\n",rate[r]); if(rpart[r]!='?') printf(" product type, parameter: %c %lf\n",rpart[r],rpar[r]); printf(" probability of reaction per time step: %lf\n",rate2[r]<0?0:rate2[r]); printf(" actual rate constant: %lf\n",calcrate(sim,i,0,r,NULL)); if(nprod[r]==2&&sim->rxn[2]&&sim->rxn[2]->nrxn[i3=prod[r][0]->ident*maxident+prod[r][1]->ident]) for(j3=0;j3rxn[2]->nrxn[i3];j3++) { r3=sim->rxn[2]->table[i3][j3]; if(sim->rxn[2]->rate2[r3]>=0) { bval=distanceVVD(prod[r][0]->pos,prod[r][1]->pos,dim); dsum=sim->mols->difc[prod[r][0]->ident]+sim->mols->difc[prod[r][1]->ident]; step=sqrt(2.0*sim->dt*dsum); pgem=1.0-numrxnrate(step,sqrt(sim->rxn[2]->rate2[r3]),-1)/numrxnrate(step,sqrt(sim->rxn[2]->rate2[r3]),bval); rev=sim->rxn[2]->nprod[r3]==1&&sim->rxn[2]->prod[r3][0]->ident==i; printf(" prob. of geminate %s reaction '%s' is %lf\n",rev?"reverse":"continuation",sim->rxn[2]->rname[r3],pgem); }}}} else if(order==2) { for(i1=0;i1 ",rname[r],name[i1],name[i2]); for(p=0;pident],p+1pos,prod[r][p]->pos,dim)>0) { printf(" product %s at ",name[prod[r][p]->ident]); printVD(prod[r][p]->pos,dim); } if(rate[r]>=0) printf(" requested rate constant: %lf\n",rate[r]); if(rpart[r]!='?') printf(" product type, parameter: %c %lf\n",rpart[r],rpar[r]); printf(" binding radius: %lf\n",rate2[r]<0?0:sqrt(rate2[r])); printf(" actual rate constant: %lf\n",calcrate(sim,i1,i2,r,NULL)); if(rate2[r]>=0) printf(" diffusion limited rate constant: %lf\n",4*PI*dsum*sqrt(rate2[table[i][j]])); if(rate2[r]>=0) printf(" activation limited rate constant: %lf\n",actrxnrate(sqrt(2.0*dsum*sim->dt),sqrt(rate2[table[i][j]]))/sim->dt); if(nprod[r]==2&&nrxn[i3=prod[r][0]->ident*maxident+prod[r][1]->ident]) for(j3=0;j3rate2[r3]>=0) { bval=distanceVVD(prod[r][0]->pos,prod[r][1]->pos,dim); dsum=sim->mols->difc[prod[r][0]->ident]+sim->mols->difc[prod[r][1]->ident]; step=sqrt(2.0*sim->dt*dsum); pgem=1.0-numrxnrate(step,sqrt(rate2[r3]),-1)/numrxnrate(step,sqrt(rate2[r3]),bval); rev=(nprod[r3]==2)&&((prod[r3][0]->ident==i1&&prod[r3][1]->ident==i2)||(prod[r3][0]->ident==i2&&prod[r3][1]->ident==i1)); printf(" prob. of geminate %s reaction '%s' is %lf\n",rev?"reverse":"continuation",sim->rxn[2]->rname[r3],pgem); }}}}} printf("\n"); return; } /* findreverserxn inputs the reaction defined by reactants i1 and/or i2 and reaction number r and looks to see if there is a reverse reaction. If both i1 and i2 are 0, then the forward reaction is zeroth order and there is no direct reverse reaction. If either i1 or i2 is 0, the forward reaction is first order, and if neither reactant is 0, the forward reaction is second order. If there is no reaction number r, an error code of Ð1 is returned. If there is a direct reverse reaction, meaning the products of the input reaction are themselves able to react to form identities i1 and i2 (or just one of them if the input reaction is first order), then the function returns 1 and the order and reaction number of the reverse reaction are pointed to by optr and rptr. If there is no direct reverse reaction, but the products of the input reaction are still able to react, the function returns 2 and optr and rptr point to the first listed continuation reaction. If the products do not react, the function returns 0 (this also includes the situation where the products of the input reaction are A and B and there is no A+B reaction, although A and/or B can undergo unimolecular reactions, as well as all situations in which there are three or more products of a forward reaction). Either or both of optr and rptr are allowed to be sent in as NULL values if the respective pieces of output information are not of interest. */ int findreverserxn(simptr sim,int i1,int i2,int r,int *optr,int *rptr) { int rev,order,o2,r2,i3,i4,i5,j2; rxnptr rxn,rxn2; if(!sim||i1<0||i2<0||i1>=sim->maxident||i2>=sim->maxident||r<0) return -1; o2=r2=rev=0; if(!i1&&!i2) order=0; else if(i1&&!i2) order=1; else if(i2&&!i1) { order=1; i1=i2; i2=0; } else order=2; rxn=sim->rxn[order]; if(!rxn||r>=rxn->total) return -1; o2=rxn->nprod[r]; if(o2==0||o2>2||!sim->rxn[o2]) return 0; rxn2=sim->rxn[o2]; i3=rxn->prod[r][0]->ident; if(o2==2) i3=sim->maxident*i3+rxn->prod[r][1]->ident; if(order>0) for(j2=0;j2nrxn[i3]&&!rev;j2++) { r2=rxn2->table[i3][j2]; if(rxn2->nprod[r2]==order) { i4=rxn2->prod[r2][0]->ident; i5=order==2?rxn2->prod[r2][1]->ident:0; rev=((i4==i1&&i5==i2)||(i4==i2&&i5==i1)); }} if(!rev&&rxn2->nrxn[i3]) { r2=0; rev=2; } if(optr) *optr=o2; if(rptr) *rptr=r2; return rev; } /* setrates is used to convert the requested reaction rates to values that are useful for the simulation routines. Values in the rate element are read and values are written to rate2. If a rate element is less than zero, it is assumed to have been unassigned and is skipped; in this case, the respective value for rate2 is not modified. For zeroth order reactions, rate2 is the expectation number of molecules that should be produced in the entire simulation volume during one time step, which is rate*dt*volume. For first order reactions, rate2 is the probability of a unimolecular reaction occuring for an individual reactant molecule during one time step, which is rate/sum*[1Ðexp(Ðsum*dt)], where sum is the sum of the defined rate values for all unimolecular reactions of the reactant. For second order reactions, rate2 is the squared binding radius of the reactants, found from bindingradius. In this case, the reverse parameter is accounted for in the reaction rate calculation if there is a direct reverse reaction and if it is appropriate (see the discussion of ÒBinding and unbinding radii,Ó above and the description for findreverserxn). This routine also sets the lists member of reaction structures. The return value of the function is Ð1 for correct operation. If errors occur, which is only possible for illegal inputs (return value of 0) or bimolecular calculations, the reaction number where the error was encountered is returned. Other than illegal inputs, the only possible errors arise from a diffusion constant of 0 for both reactants, or a directly reversible reaction that has an undeclared reversible type. */ int setrates(simptr sim,int order) { int dim,nident,maxident; int r,i,i1,i2,j,d,rev,o2,r2,ll1,ll2; double dsum,*difc,dt,rate3,rpar,flt; double *rate,*rate2; rxnptr rxn; char rpart; if(!sim||order<0||order>2) return 0; rxn=sim->rxn[order]; if(!rxn) return -1; dim=sim->dim; maxident=sim->maxident; nident=sim->nident; difc=sim->mols->difc; dt=sim->dt; rate=rxn->rate; rate2=rxn->rate2; if(order==0) { // order 0 flt=1.0; for(d=0;dwlist[2*d+1]->pos-sim->wlist[2*d]->pos; for(r=0;rtotal;r++) if(rate2[r]<0&&rate[r]>=0) rate2[r]=rate[r]*dt*flt; } else if(order==1) { // order 1 for(i=0;inrxn[i];j++) { r=rxn->table[i][j]; flt+=rate[r]>=0?rate[r]:0; } for(j=0;jnrxn[i];j++) { r=rxn->table[i][j]; if(rate2[r]>=0); else if(rate[r]==0) rate2[r]=0; else if(rate[r]>0) rate2[r]=rate[r]/flt*(1.0-exp(-dt*flt)); }}} else if(order==2) { // order 2 for(i1=0;i1nrxn[i];j++) { r=rxn->table[i][j]; if(i1==i2) rate3=2*rate[r]; else rate3=rate[r]; dsum=difc[i1]+difc[i2]; rev=findreverserxn(sim,i1,i2,r,&o2,&r2); if(rev==1) { rpart=sim->rxn[o2]->rpart[r2]; rpar=sim->rxn[o2]->rpar[r2]; } else { rpart=' '; rpar=0; } if(rate[r]<0); else if(rate2[r]>=0); else if(rate3==0) rate2[r]=0; else if(dsum<=0) return r; else if(rev!=1) rate2[r]=bindingradius(rate3,dt,dsum,-1,0); else if(strchr("qysofi",rpart)) rate2[r]=bindingradius(rate3,dt,dsum,-1,0); else if(rpart=='b') rate2[r]=bindingradius(rate3,dt,dsum,rpar,0); else if(rpart=='r') rate2[r]=bindingradius(rate3,dt,dsum,rpar,1); else if(rpart=='p') rate2[r]=bindingradius(rate3*(1.0-rpar),dt,dsum,-1,0); else if(rpart=='x') { rate2[r]=bindingradius(rate3,dt,dsum,0,0); flt=unbindingradius(rpar,dt,dsum,rate2[r]); if(flt>0) rate2[r]=bindingradius(rate3*(1-rpar),dt,dsum,-1,0); } else return r; if(rate2[r]>=0) rate2[r]*=rate2[r]; }}} if(order==1) { // set rxn->lists information for(i=1;inrxn[i];r++) if(rate2[r]>0) rxn->lists|=(sim->mols->difc[i]==0)+1; } else if(order==2) { for(i1=1;i1nrxn[maxident*i1+i2];r++) { ll1=sim->mols->difc[i1]==0; ll2=sim->mols->difc[i2]==0; if(rate2[r]>0) { if(ll1==0&&ll2==0) rxn->lists|=1; else if((ll1==0&&ll2==1)||(ll1==1&&ll2==0)) rxn->lists|=2; }}} return -1; } /* setproducts is used to set the initial separations between reaction products for all reactions of order order, based on the corresponding rpart character and rpar parameter. This is done by calculating the proper separation and then setting the first value of the template molecule pos vector to this separation, for all appropriate templates. If rpart is either ÔoÕ or ÔfÕ, denoting either randomly oriented offset or fixed orientation offset, then it is assumed that the template molecule position has already been set up; it is not modified again by this routine. Otherwise, it is assumed that the template molecule position vectors have all values equal to 0 initially. If there are illegal inputs, 0 is returned by this routine. If an error occurs, the reaction number where the error was encountered is returned and a message is returned in the string erstr, which should have been allocated to size STRCHAR; if all assignments work correctly, Ð1 is returned and the string is unchanged; and if a reversible reaction was undeclared, Ð1 is returned and a warning is returned in the string, although the program does not need to terminate. Possible errors include trying to set product positions for reactions without products, setting relative positions for products in which reactions have only one product, are irreversible, or have multiple reverse reactions, or trying to get a geminate binding probability that is unachievably high due to too long a time step. See the discussion in the section called ÒBinding and unbinding radiiÓ for more details. */ int setproducts(simptr sim,int order,char *erstr) { int maxident,nident,i1,i2,i,j,r,r2,rev,o2; rxnptr rxn,rxn2; double a,rpar,*difc,dc1,dc2,dsum; moleculeptr mptr1,mptr2; char rpart; if(!sim||order<0||order>2||!erstr) return 0; maxident=sim->maxident; nident=sim->nident; rxn=sim->rxn[order]; difc=sim->mols->difc; if(!rxn) return -1; for(i1=0;i1==0||(order>0&&i1total)||(order>0&&jnrxn[i]);j++) { if(order==0) r=j; else r=rxn->table[i][j]; rpar=rxn->rpar[r]; rpart=rxn->rpart[r]; o2=rxn->nprod[r]; if(o2==0&&strchr("pxrbqysof",rpart)) { strncpy(erstr,"Illegal product parameter because reaction has no products",STRCHAR); return r; } else if(o2==1) { mptr1=rxn->prod[r][0]; if(strchr("pxrqys",rpart)) { strncpy(erstr,"Illegal product parameter because reaction only has one product",STRCHAR); return r; } else if(rpart=='b') mptr1->pos[0]=rpar; } else if(o2==2) { mptr1=rxn->prod[r][0]; mptr2=rxn->prod[r][1]; dc1=difc[mptr1->ident]; dc2=difc[mptr2->ident]; dsum=dc1+dc2; rev=findreverserxn(sim,i1,i2,r,&o2,&r2); if(rev>0) rxn2=sim->rxn[o2]; if(rev>0&&rxn2->rate2[r2]>=0) a=sqrt(rxn2->rate2[r2]); else a=-1; if(rev<=0) { if(strchr("pxrqys",rpart)) { strncpy(erstr,"Illegal product parameter because products don't react",STRCHAR); return r; } else if(strchr("ofi?",rpart)); else if(dsum<=0) { strncpy(erstr,"Cannot set unbinding distance because sum of product diffusion constants is 0",STRCHAR); return r; } else { mptr1->pos[0]=rpar*dc1/dsum; mptr2->pos[0]=-rpar*dc2/dsum; }} else { if(strchr("ofi",rpart)); else if(rpart=='?') strncpy(erstr,"WARNING: unbinding radius was set to 0 because of an undefined product parameter",STRCHAR); else if(dsum<=0) { strncpy(erstr,"Cannot set unbinding distance because sum of product diffusion constants is 0",STRCHAR); return r; } else if(rpart=='b') { mptr1->pos[0]=rpar*dc1/dsum; mptr2->pos[0]=-rpar*dc2/dsum; } else if(rxn2->rate2[r2]<0) strncpy(erstr,"WARNING: Binding radius of reaction products is undefined",STRCHAR); else if(strchr("rs",rpart)) { mptr1->pos[0]=rpar*a*dc1/dsum; mptr2->pos[0]=-rpar*a*dc2/dsum; } else if(strchr("pq",rpart)) { rpar=unbindingradius(rpar,sim->dt,dsum,a); if(rpar==-2) { strncpy(erstr,"Cannot create an unbinding radius due to illegal input values",STRCHAR); return r; } else if(rpar<0) { strncpy(erstr,"Maximum possible geminate binding probability is ",STRCHAR); sprintf(erstr,"%lf",-rpar); return r; } else { mptr1->pos[0]=rpar*dc1/dsum; mptr2->pos[0]=-rpar*dc2/dsum; }} else if(strchr("xy",rpart)) { rpar=unbindingradius(rpar,sim->dt,dsum,a); if(rpar==-2) strncpy(erstr,"WARNING: Possible illegal input values",STRCHAR); else if(rpar>0) { mptr1->pos[0]=rpar*dc1/dsum; mptr2->pos[0]=-rpar*dc2/dsum; }}}}}} return -1; } /* calcrate calculates the macroscopic rate constant using the microscopic parameters that have been calculated or that were initially assigned. All going well, these results should exactly match those that were requested initially, although this routine is useful as a check, and for situations where the microscopic values were input rather than the mass action rate constants. For bimolecular reactions that are reversible, the routine calculates rates with accounting for reversibility if the product parameter is p, x, r, b, or ?, and not otherwise. A value of Ð1 is returned if input parameters are illegal and a value of 0 is returned if the rate2 value for the indicated reaction is undefined (<0). If reversibility is accounted for and pgemptr is not input as NULL, *pgemptr is set to the probability of geminate recombination of the reactants; otherwise its value is not changed. */ double calcrate(simptr sim,int i1,int i2,int r,double *pgemptr) { double rate3; int order,dim,maxident,nident; int i,j,d,r2,rev,o2; double flt,flt2,step,*difc,dt,a,bval; double *rate2; rxnptr rxn; moleculeptr **prod; if(!sim) return -1; if(!i1&&!i2) order=0; else if(i1&&!i2) order=1; else if(i2&&!i1) { order=1; i1=i2; i2=0; } else order=2; rxn=sim->rxn[order]; if(!rxn||r<0||r>=rxn->total) return -1; dim=sim->dim; maxident=sim->maxident; nident=sim->nident; difc=sim->mols->difc; dt=sim->dt; rate2=rxn->rate2; if(rate2[r]<0) rate3=0; else if(order==0) { // order 0 flt=1.0; for(d=0;dwlist[2*d+1]->pos-sim->wlist[2*d]->pos; rate3=rate2[r]/dt/flt; } else if(order==1) { // order 1 flt=0; for(j=0;jnrxn[i1];j++) { r2=rxn->table[i1][j]; flt+=rate2[r2]>=0?rate2[r2]:0; } flt2=-log(1.0-flt)/dt; rate3=rate2[r]>0?rate2[r]/flt*flt2:0; } else if(order==2) { // order 2 i=maxident*i1+i2; step=sqrt(2.0*(difc[i1]+difc[i2])*dt); a=sqrt(rate2[r]); rev=findreverserxn(sim,i1,i2,r,&o2,&r2); if(rev!=1||strchr("qysofi",rxn->rpart[r])) rate3=numrxnrate(step,a,-1); else { prod=sim->rxn[o2]->prod; bval=distanceVVD(prod[r2][0]->pos,prod[r2][1]->pos,dim); rate3=numrxnrate(step,a,bval); if(pgemptr) *pgemptr=1.0-numrxnrate(step,a,-1)/rate3; } rate3/=dt; if(i1==i2) rate3/=2.0; } else rate3=0; return rate3; } /* doreact executes a reaction that has already been determined to have happened. rxn is the reaction, r is the reaction number and mptr1 and mptr2 are the reactants, where mptr2 is ignored for unimolecular reactions, and both are ignored for zeroth order reactions. Reactants are killed, but left in the live lists. Any products are created on the dead list, for transfer to the live list by the molsort routine. Molecules that are created are put at the reaction position, which is the average position of the reactants weighted by the inverse of their diffusion constants, plus an offset from the product definition. The cluster of products is typically rotated to a random orientation. If the displacement was set to all 0Õs (recommended for non-reacting products), the routine is fairly fast, putting all products at the reaction position. If the rpart character is ÔfÕ, the orientation is fixed and there is no rotation. Otherwise, a non-zero displacement results in the choosing of random angles and vector rotations. If the system has more than three dimensions, only the first three are randomly oriented, while higher dimensions just add the displacement to the reaction position. The function returns 0 for successful operation and 1 if more molecules are required than were initially allocated. */ int doreact(rxnptr rxn,int r,moleculeptr mptr1,moleculeptr mptr2,simptr sim) { int p,d,nprod,dim; int rot,calc; double dc1,dc2; molssptr mols; moleculeptr mptr,mptrt; boxptr bptr; wallptr *wlist; double *v1,*v2,*m3; mols=sim->mols; dim=sim->dim; v1=sim->v1; v2=sim->v2; m3=sim->m3; if(rxn->order==0) { // kill reactants and get reaction position in v2 wlist=sim->wlist; for(d=0;dpos,wlist[2*d+1]->pos); bptr=pos2box(sim,v2); } else if(rxn->order==1) { for(d=0;dpos[d]; mptr1->ident=0; bptr=mptr1->box; } else if(rxn->order==2) { dc1=mols->difc[mptr1->ident]; dc2=mols->difc[mptr2->ident]; mptr1->ident=0; mptr2->ident=0; if(dc1==0&&dc2==0) dc1=dc2=1; sumVD(dc2/(dc1+dc2),mptr1->pos,dc1/(dc1+dc2),mptr2->pos,v2,dim); bptr=mptr1->box; } else bptr=NULL; nprod=rxn->nprod[r]; // place products calc=0; for(p=0;ptopd) return 1; mptr=mols->dead[--mols->topd]; mptrt=rxn->prod[r][p]; mptr->ident=mptrt->ident; mptr->box=bptr; for(d=0;dposx[d]=v2[d]; for(d=0;dwrap[d]=0; rot=0; for(d=0;dpos[d]!=0; if(!rot) for(d=0;dpos[d]=v2[d]; else if(rxn->rpart[r]=='f') sumVD(1,v2,1,mptrt->pos,mptr->pos,dim); else if(dim==1) { if(!calc) {m3[0]=signrand();calc=1;} mptr->pos[0]=v2[0]+m3[0]*mptrt->pos[0]; } else if(dim==2) { if(!calc) {DirCosM2D(m3,unirand30(0,2*PI));calc=1;} dotMVD(m3,mptrt->pos,v1,2,2); sumVD(1,v2,1,v1,mptr->pos,2); } else if(dim==3) { if(!calc) {DirCosMD(m3,thetarand(),unirand30(0,2*PI),unirand30(0,2*PI));calc=1;} dotMVD(m3,mptrt->pos,v1,3,3); sumVD(1,v2,1,v1,mptr->pos,3);} else { if(!calc) {DirCosMD(m3,thetarand(),unirand30(0,2*PI),unirand30(0,2*PI));calc=1;} dotMVD(m3,mptrt->pos,v1,3,3); for(d=3;dpos[d]; sumVD(1,v2,1,v1,mptr->pos,dim); }} return 0; } /* zeroreact figures out how many molecules to create for each zeroth order reaction and then tells doreact to create them. It returns the number of molecules created, or Ð1 if not enough molecules were allocated initially. */ int zeroreact(simptr sim) { int r,rctr,nmol; rxnptr rxn; rxn=sim->rxn[0]; if(!rxn) return 0; rctr=0; for(r=0;rtotal;r++) { nmol=poisrandD(rxn->rate2[r]); rctr+=nmol; for(;nmol>0;nmol--) if(doreact(rxn,r,NULL,NULL,sim)) return -1; } return rctr; } /* unireact determines whether unimolecular reactions occured, considering both live lists. Reactions that do occur are sent to doreact to process them. The function returns the number of reactions that occured during that time step, or Ð1 if not enough molecules were allocated initially. */ int unireact(simptr sim) { rxnptr rxn; moleculeptr *mlist; int rctr,*nrxn,**table; double *rate2; int i,j,r,m,nmol,ll; rxn=sim->rxn[1]; if(!rxn) return 0; rctr=0; nrxn=rxn->nrxn; table=rxn->table; rate2=rxn->rate2; for(ll=0;ll<2;ll++) if(rxn->lists&(ll+1)) { mlist=sim->mols->live[ll]; nmol=sim->mols->nl[ll]; for(m=0;mident; for(j=0;jrxn[2]; if(!rxn) return 0; dim=sim->dim; live=sim->mols->live; maxident=sim->maxident; rctr=0; nrxn=rxn->nrxn; table=rxn->table; rate2=rxn->rate2; nmol=sim->mols->nl; if(!nmol[0]) return 0; if(!neigh) { // same box for(ll=0;ll<2;ll++) if(rxn->lists&(ll+1)) for(m1=0;m1box; mlist2=bptr->mol[0]; nmol2=bptr->nmol[0]; for(m2=0;m2ident*maxident+mptr2->ident; if(nrxn[i]) { for(dist2=d=0;dpos[d]-mptr2->pos[d])*(mptr1->pos[d]-mptr2->pos[d]); for(j=0;jlists&(ll+1)) for(m1=0;m1box; bmax=ll?bptr->nneigh:bptr->midneigh; for(b2=0;b2neigh[b2]->mol[0]; nmol2=bptr->neigh[b2]->nmol[0]; if(bptr->wpneigh&&bptr->wpneigh[b2]) { wpcode=bptr->wpneigh[b2]; for(m2=0;m2ident*maxident+mptr2->ident; if(nrxn[i]) { for(dist2=d=0;d>2*d&3)==0) dist2+=(mptr1->pos[d]-mptr2->pos[d])*(mptr1->pos[d]-mptr2->pos[d]); else if((wpcode>>2*d&3)==1) { pos2=sim->wlist[2*d+1]->pos-sim->wlist[2*d]->pos; dist2+=(mptr1->pos[d]-mptr2->pos[d]+pos2)*(mptr1->pos[d]-mptr2->pos[d]+pos2); } else { pos2=sim->wlist[2*d+1]->pos-sim->wlist[2*d]->pos; dist2+=(mptr1->pos[d]-mptr2->pos[d]-pos2)*(mptr1->pos[d]-mptr2->pos[d]-pos2); } for(j=0;jmols->difc[mptr1->ident]; dc2=sim->mols->difc[mptr2->ident]; if(dc1pos[d]=mptr1->pos[d]; else for(d=0;dpos[d]=mptr2->pos[d]; if(doreact(rxn,table[i][j],mptr1,mptr2,sim)) return -1; }}}} else for(m2=0;m2ident*maxident+mptr2->ident; if(nrxn[i]) { for(dist2=d=0;dpos[d]-mptr2->pos[d])*(mptr1->pos[d]-mptr2->pos[d]); for(j=0;j1) { ps=3; npts=3; } else if(pshape=='h'&&dim>1) { ps=4; npts=3; } else { ps=-1; npts=0; } if(nptsptr) *nptsptr=npts; return ps; } /* Allocates maxpanels of shape pshape for the surface srf. The srf element of the panels are set to srf. In srf, the correct maxpanel entry is set to maxpanel, the npanel entry is set to 0, and the proper list of panels are allocated and cleared. All points are set to all zeros. The function returns 1 for success and 0 for failure to allocate memory. */ int panelsalloc(surfaceptr srf,int dim,int maxpanel,char pshape) { panelptr *pnls,pnl; int p,npts,pt,d,ps; ps=pshape2ps(pshape,dim,&npts); if(ps<0) return 0; pnls=(panelptr*) calloc(maxpanel,sizeof(panelptr)); if(!pnls) return 0; for(p=0;ppshape=pshape; pnl->srf=srf; pnl->npts=npts; pnl->point=NULL; CHECK(pnl->point=(double**) calloc(npts,sizeof(double*))); for(pt=0;ptpoint[pt]=NULL; for(pt=0;ptpoint[pt]=(double*) calloc(dim,sizeof(double))); for(d=0;dpoint[pt][d]=0; } pnl->front[0]=pnl->front[1]=pnl->front[2]=0; } if(srf) { srf->maxpanel[ps]=maxpanel; srf->npanel[ps]=0; srf->panels[ps]=pnls; } return 1; failure: for(p=0;pnpts&&pnl->point) { for(pt=0;ptnpts;pt++) if(pnl->point[pt]) free(pnl->point[pt]); free(pnl->point); } free(pnl); return; } /* Allocates a surface structure, and sets all elements to initial values. maxident is the number of molecular identities, which is used for allocating faction and baction. dim is the system dimensionality, as always. Colors are set to all 0Õs (black), but with alpha values of 1 (opaque); polygon modes are set to ÔfÕ for face; edgepoints is set to 1; faction and baction are set to ÔrÕ for reflective. No panels are allocated. This is called by surfacessalloc and so should not need to be called externally. */ surfaceptr surfacealloc(int maxident,int dim) { surfaceptr srf; int i,ps; srf=(surfaceptr) malloc(sizeof(struct surfacestruct)); if(!srf) return NULL; srf->faction=NULL; srf->baction=NULL; srf->fcolor[0]=srf->fcolor[1]=srf->fcolor[2]=0; srf->bcolor[0]=srf->bcolor[1]=srf->bcolor[2]=0; srf->fcolor[3]=srf->bcolor[3]=1; srf->edgepts=1; srf->fpolymode=(dim==3?'f':'e'); srf->bpolymode=(dim==3?'f':'e'); for(ps=0;psmaxpanel[ps]=0; for(ps=0;psnpanel[ps]=0; for(ps=0;pspanels[ps]=NULL; srf->faction=(char*) calloc(maxident,sizeof(char)); if(!srf->faction) {surfacefree(srf);return NULL;} for(i=0;ifaction[i]='r'; srf->baction=(char*) calloc(maxident,sizeof(char)); if(!srf->baction) {surfacefree(srf);return NULL;} for(i=0;ibaction[i]='r'; return srf; } /* Frees a surface, including all substructures and panels in it. This is called by surfacessfree and so should not need to be called externally. */ void surfacefree(surfaceptr srf) { int p,ps; if(!srf) return; for(ps=0;psmaxpanel[ps];p++) panelfree(srf->panels[ps][p]); free(srf->panels[ps]); } free(srf->faction); free(srf->baction); free(srf); return; } /* Allocates a surface superstructure for maxsurface surfaces, as well as all of the surfaces. Each name is allocated to an empty string of STRCHAR (256) characters. */ surfacessptr surfacessalloc(int maxsurface,int maxident,int dim) { surfacessptr srfss; int s; srfss=(surfacessptr) malloc(sizeof(struct surfacesuperstruct)); if(!srfss) return NULL; srfss->maxsrf=maxsurface; srfss->nsrf=0; srfss->snames=NULL; srfss->srflist=NULL; CHECK(srfss->snames=(char**) calloc(maxsurface,sizeof(char*))); for(s=0;ssnames[s]=NULL; for(s=0;ssnames[s]=EmptyString()); CHECK(srfss->srflist=(surfaceptr*) calloc(maxsurface,sizeof(surfaceptr))); for(s=0;ssrflist[s]=NULL; for(s=0;ssrflist[s]=surfacealloc(maxident,dim)); return srfss; failure: surfacessfree(srfss); return NULL; } /* Frees a surface superstructure pointed to by srfss, and all contents in it, including all of the surfaces and all of their panels. */ void surfacessfree(surfacessptr srfss) { int s; if(!srfss) return; if(srfss->srflist) { for(s=0;smaxsrf;s++) surfacefree(srfss->srflist[s]); free(srfss->srflist); } if(srfss->snames) { for(s=0;smaxsrf;s++) free(srfss->snames[s]); free(srfss->snames); } free(srfss); return; } /* loadsurface loads a surface from an already opened disk file pointed to with fptr. lctrptr is a pointer to the line counter, which is updated each time a line is read. If successful, it returns 0 and the surface is added to the surface superstructure in sim, which should have been already allocated. Otherwise it returns the updated line counter along with an error message. If a surface with the same name (entered by the user) already exists, this function can add more panels to it. It can also allocate and set up a new surface. If this runs successfully, the complete surface structure is set up, with the exception of box issues. If the routine fails, any new surface structure is freed. */ int loadsurface(simptr sim,FILE *fptr,int *lctrptr,char *erstr) { surfacessptr srfss; surfaceptr srf; char line[STRCHAR],word[STRCHAR],*line2; char nm[STRCHAR],ch,pshape; int got[6],itct,dim,incomment; int s,i1,pdim,pt,ps,i; panelptr pnl; double fltv1[9],**point,*front; setstdZV(got,6,0); incomment=0; srfss=sim->srfss; CHECKS(srfss,"PROGRAM BUG: surface superstructure not allocated in loadsurface"); srf=NULL; dim=sim->dim; while(fgets(line,STRCHAR,fptr)) { (*lctrptr)++; if(strchr(line,'#')) *strchr(line,'#')='\0'; itct=sscanf(line,"%s",word); if(itct) line2=strnword(line,2); else line2=NULL; if(itct<=0); // blank line else if(!strcmp(word,"/*")) { // start comment with /* incomment=1; } else if(!strcmp(word,"*/")) { // end comment with */ CHECKS(incomment,"*/ without a corresponding /*"); incomment=0; } else if(incomment); else if(!strcmp(word,"end_surface")) { // end_surface CHECKS(!line2,"unexpected text following end_surface"); return 0; } else if(!line2) { // just word CHECKS(0,"unknown word or missing parameter"); } else if(!strcmp(word,"name")) { // name, got[0] CHECKS(!got[0],"name can only be entered once in each surface block"); got[0]++; itct=sscanf(line2,"%s",nm); CHECKS(itct==1,"error reading surface name"); s=stringfind(srfss->snames,srfss->nsrf,nm); if(s<0) { CHECKS(srfss->nsrfmaxsrf,"more surfaces are being defined than were allocated"); s=srfss->nsrf++; strncpy(srfss->snames[s],nm,STRCHAR); srf=srfss->srflist[s]; } else srf=srfss->srflist[s]; } else if(!strcmp(word,"action_front")) { // action_front CHECKS(got[0],"need to enter surface name before action_front"); itct=sscanf(line2,"%s %c",nm,&ch); CHECKS(itct==2,"action_front format: molecule_name action"); if(!strcmp(nm,"all")) i=-2; else i=stringfind(sim->name,sim->nident,nm); CHECKS(i!=-1,"in action_front, molecule name not recognized"); CHECKS(ch=='r'||ch=='a'||ch=='t'||ch=='p',"action_front actions allowed are 'r', 'a', 't', or 'p'"); if(i==-2) for(i=0;inident;i++) srf->faction[i]=ch; else srf->faction[i]=ch; CHECKS(!strnword(line2,3),"unexpected text following action_front"); } else if(!strcmp(word,"action_back")) { // action_back CHECKS(got[0],"need to enter surface name before action_back"); itct=sscanf(line2,"%s %c",nm,&ch); CHECKS(itct==2,"action_back format: molecule_name action"); if(!strcmp(nm,"all")) i=-2; else i=stringfind(sim->name,sim->nident,nm); CHECKS(i!=-1,"in action_back, molecule name not recognized"); CHECKS(ch=='r'||ch=='a'||ch=='t'||ch=='p',"action_back actions allowed are 'r', 'a', 't', or 'p'"); if(i==-2) for(i=0;inident;i++) srf->baction[i]=ch; else srf->baction[i]=ch; CHECKS(!strnword(line2,3),"unexpected text following action_back"); } else if(!strcmp(word,"action_both")) { // action_both CHECKS(got[0],"need to enter surface name before action_both"); itct=sscanf(line2,"%s %c",nm,&ch); CHECKS(itct==2,"action_both format: molecule_name action"); if(!strcmp(nm,"all")) i=-2; else i=stringfind(sim->name,sim->nident,nm); CHECKS(i!=-1,"in action_both, molecule name not recognized"); CHECKS(ch=='r'||ch=='a'||ch=='t'||ch=='p',"action_both actions allowed are 'r', 'a', 't', or 'p'"); if(i==-2) for(i=0;inident;i++) srf->faction[i]=srf->baction[i]=ch; else srf->faction[i]=srf->baction[i]=ch; CHECKS(!strnword(line2,3),"unexpected text following action_both"); } else if(!strcmp(word,"color_front")) { // color_front, got[1] CHECKS(got[0],"need to enter surface name before color_front"); CHECKS(!got[1],"color_front can only be entered once for each surface"); got[1]++; itct=sscanf(line2,"%lf %lf %lf",&(srf->fcolor[0]),&(srf->fcolor[1]),&(srf->fcolor[2])); CHECKS(itct==3,"surface color format: red green blue [alpha]"); CHECKS(srf->fcolor[0]>=0&&srf->fcolor[0]<=1,"color values need to between 0 and 1"); CHECKS(srf->fcolor[1]>=0&&srf->fcolor[1]<=1,"color values need to between 0 and 1"); CHECKS(srf->fcolor[2]>=0&&srf->fcolor[2]<=1,"color values need to between 0 and 1"); if((line2=strnword(line2,4))!=NULL) { itct=sscanf(line2,"%lf",&(srf->fcolor[3])); CHECKS(srf->fcolor[3]>=0&&srf->fcolor[3]<=1,"alpha value needs to be between 0 and 1"); line2=strnword(line2,2); } CHECKS(!line2,"unexpected text following color_front"); } else if(!strcmp(word,"color_back")) { // color_back, got[2] CHECKS(got[0],"need to enter surface name before color_back"); CHECKS(!got[2],"color_back can only be entered once for each surface"); got[2]++; itct=sscanf(line2,"%lf %lf %lf",&(srf->bcolor[0]),&(srf->bcolor[1]),&(srf->bcolor[2])); CHECKS(itct==3,"surface color format: red green blue"); CHECKS(srf->bcolor[0]>=0&&srf->bcolor[0]<=1,"color values need to between 0 and 1"); CHECKS(srf->bcolor[1]>=0&&srf->bcolor[1]<=1,"color values need to between 0 and 1"); CHECKS(srf->bcolor[2]>=0&&srf->bcolor[2]<=1,"color values need to between 0 and 1"); if((line2=strnword(line2,4))!=NULL) { itct=sscanf(line2,"%lf",&(srf->bcolor[3])); CHECKS(srf->bcolor[3]>=0&&srf->bcolor[3]<=1,"alpha value needs to be between 0 and 1"); line2=strnword(line2,2); } CHECKS(!line2,"unexpected text following color_back"); } else if(!strcmp(word,"color")||!strcmp(word,"color_both")) { // color, color_both, got[1,2] CHECKS(got[0],"need to enter surface name before color"); CHECKS(!got[1]&&!got[2],"colors can only be entered once for each side, for each surface"); got[1]++; got[2]++; itct=sscanf(line2,"%lf %lf %lf",&(srf->fcolor[0]),&(srf->fcolor[1]),&(srf->fcolor[2])); CHECKS(itct==3,"surface color format: red green blue [alpha]"); CHECKS(srf->fcolor[0]>=0&&srf->fcolor[0]<=1,"color values need to between 0 and 1"); CHECKS(srf->fcolor[1]>=0&&srf->fcolor[1]<=1,"color values need to between 0 and 1"); CHECKS(srf->fcolor[2]>=0&&srf->fcolor[2]<=1,"color values need to between 0 and 1"); if((line2=strnword(line2,4))!=NULL) { itct=sscanf(line2,"%lf",&(srf->fcolor[3])); CHECKS(srf->fcolor[3]>=0&&srf->fcolor[3]<=1,"alpha value needs to be between 0 and 1"); line2=strnword(line2,2); } srf->bcolor[0]=srf->fcolor[0]; srf->bcolor[1]=srf->fcolor[1]; srf->bcolor[2]=srf->fcolor[2]; srf->bcolor[3]=srf->fcolor[3]; CHECKS(!line2,"unexpected text following color"); } else if(!strcmp(word,"thickness")) { // thickness, got[3] CHECKS(got[0],"need to enter surface name before thickness"); CHECKS(!got[3],"thickness can only be entered once for each surface"); got[3]++; itct=sscanf(line2,"%lf",&(srf->edgepts)); CHECKS(itct==1,"thickness value is missing"); CHECKS(srf->edgepts>=0,"thickness value needs to be at least 0"); CHECKS(!strnword(line2,2),"unexpected text following thickness"); } else if(!strcmp(word,"polygon_front")) { // polygon_front, got[4] CHECKS(got[0],"need to enter surface name before polygon_front"); CHECKS(!got[4],"polygon_front can only be entered once for each surface"); got[4]++; itct=sscanf(line2,"%c",&ch); CHECKS(itct==1,"polygon_front character is missing"); CHECKS(ch=='f'||ch=='e'||ch=='g'||ch=='v',"polygon_front needs to be 'f', 'e', 'g', or 'v'"); srf->fpolymode=ch; CHECKS(!strnword(line2,2),"unexpected text following polygon_front"); } else if(!strcmp(word,"polygon_back")) { // polygon_back, got[5] CHECKS(got[0],"need to enter surface name before polygon_back"); CHECKS(!got[5],"polygon_back can only be entered once for each surface"); got[5]++; itct=sscanf(line2,"%c",&ch); CHECKS(itct==1,"polygon_back character is missing"); CHECKS(ch=='f'||ch=='e'||ch=='g'||ch=='v',"polygon_back needs to be 'f', 'e', 'g', or 'v'"); srf->bpolymode=ch; CHECKS(!strnword(line2,2),"unexpected text following polygon_back"); } else if(!strcmp(word,"polygon_both")) { // polygon_both, got[4,5] CHECKS(got[0],"need to enter surface name before polygon_both"); CHECKS(!got[4]&&!got[5],"polygon drawing can only be entered once for each side, for each surface"); got[4]++; got[5]++; itct=sscanf(line2,"%c",&ch); CHECKS(itct==1,"polygon_both character is missing"); CHECKS(ch=='f'||ch=='e'||ch=='g'||ch=='v',"polygon_both needs to be 'f', 'e', 'g', or 'v'"); srf->fpolymode=ch; srf->bpolymode=ch; CHECKS(!strnword(line2,2),"unexpected text following polygon_both"); } else if(!strcmp(word,"max_panels")) { // max_panels CHECKS(got[0],"need to enter surface name before max_panels"); itct=sscanf(line2,"%c %i",&ch,&i1); CHECKS(itct==2,"max_panels format: shape number"); ps=pshape2ps(ch,dim,NULL); CHECKS(ps>=0,"the panel shape needs to be 'r', 't', 's', 'c', or 'h' ('c' and 'h' can't be used in 1-D)"); CHECKS(srf->maxpanel[ps]==0,"max_panels can only be entered once per surface, for each shape"); CHECKS(i1>=0,"max_panels number needs to be at least 0"); if(i1>0) CHECKS(panelsalloc(srf,dim,i1,ch),"memory error at max_panels statement"); CHECKS(!strnword(line2,3),"unexpected text following max_panels"); } else if(!strcmp(word,"panel")) { // panel CHECKS(got[0],"need to enter surface name before panel"); itct=sscanf(line2,"%c",&pshape); CHECKS(itct==1,"panel shape needs to be entered"); ps=pshape2ps(pshape,dim,NULL); CHECKS(ps>=0,"the panel shape needs to be 'r', 't', 's', 'c', or 'h' ('c' and 'h' can't be used in 1-D)"); CHECKS(srf->maxpanel[ps]>0,"need to enter max_panels before panel for this shape"); CHECKS(srf->npanel[ps]maxpanel[ps],"more panels are specified for this shape than were allocated"); pnl=srf->panels[ps][srf->npanel[ps]++]; point=pnl->point; front=pnl->front; line2=strnword(line2,2); CHECKS(line2,"panel data missing"); if(pshape=='r') { // panel r ... itct=sscanf(line2,"%c",&ch); CHECKS(itct==1,"error reading rectangle panel direction (should be + or -)"); CHECKS(ch=='+'||ch=='-',"error reading rectangle panel direction (should be + or -)"); if(ch=='+') line2=strchr(line2,'+')+1; else line2=strchr(line2,'-')+1; itct=sscanf(line2,"%i",&pdim); CHECKS(itct==1,"error reading rectangle panel dimension"); CHECKS(pdim>=0&&pdim0,"drawing slices for sphere panel needs to be positive"); front[0]=(double)sign(fltv1[2]); } else if(dim==3) { point[0][0]=fltv1[0]; point[0][1]=fltv1[1]; point[0][2]=fltv1[2]; point[1][0]=fabs(fltv1[3]); point[1][1]=fltv1[4]; CHECKS(point[1][1]>0,"drawing slices for sphere panel needs to be positive"); point[1][2]=fltv1[5]; CHECKS(point[1][2]>0,"drawing stacks for sphere panel needs to be positive"); front[0]=(double)sign(fltv1[3]); }} else if(pshape=='c') { // panel c ... itct=strreadnd(line2,4*dim-3,fltv1,NULL); CHECKS(itct==4*dim-3,"cylinder panel is missing values"); if(dim==1) { CHECKS(0,"cylinder panels cannot be used with a 1 dimensional system"); } else if(dim==2) { point[0][0]=fltv1[0]; point[0][1]=fltv1[1]; point[1][0]=fltv1[2]; point[1][1]=fltv1[3]; CHECKS(!(point[0][0]==point[1][0]&&point[0][1]==point[1][1]),"cylinder ends need to be at different locations"); point[2][0]=fabs(fltv1[4]); Geo_LineNormal(point[0],point[1],front); front[2]=(double)sign(fltv1[4]); } else if(dim==3) { point[0][0]=fltv1[0]; point[0][1]=fltv1[1]; point[0][2]=fltv1[2]; point[1][0]=fltv1[3]; point[1][1]=fltv1[4]; point[1][2]=fltv1[5]; CHECKS(!(point[0][0]==point[1][0]&&point[0][1]==point[1][1]&&point[0][2]==point[1][2]),"cylinder ends need to be at different locations"); point[2][0]=fabs(fltv1[6]); point[2][1]=fltv1[7]; CHECKS(point[2][1]>0,"drawing slices for cylinder panel needs to be positive"); point[2][2]=fltv1[8]; CHECKS(point[2][2]>0,"drawing stacks for cylinder panel needs to be positive"); front[0]=(double)sign(fltv1[6]); }} else if(pshape=='h') { // panel h ... itct=strreadnd(line2,3*dim,fltv1,NULL); CHECKS(itct==3*dim,"hemisphere panel is missing values"); if(dim==1) { CHECKS(0,"hemisphere panels cannot be used with a 1 dimensional system"); } else if(dim==2) { point[0][0]=fltv1[0]; point[0][1]=fltv1[1]; point[1][0]=fabs(fltv1[2]); point[1][1]=fltv1[5]; CHECKS(point[1][1]>0,"drawing slices for hemisphere panel needs to be positive"); point[2][0]=fltv1[3]; point[2][1]=fltv1[4]; CHECKS(normalizeVD(point[2],2)>0,"outward pointing vector cannot be 0 length"); front[0]=(double)sign(fltv1[2]); } else if(dim==3) { point[0][0]=fltv1[0]; point[0][1]=fltv1[1]; point[0][2]=fltv1[2]; point[1][0]=fabs(fltv1[3]); point[1][1]=fltv1[7]; CHECKS(point[1][1]>0,"drawing slices for hemisphere panel needs to be positive"); point[1][2]=fltv1[8]; CHECKS(point[1][2]>0,"drawing stacks for hemisphere panel needs to be positive"); point[2][0]=fltv1[4]; point[2][1]=fltv1[5]; point[2][2]=fltv1[6]; CHECKS(normalizeVD(point[2],3)>0,"outward pointing vector cannot be 0 length"); front[0]=(double)sign(fltv1[3]); }}} else { // unknown word CHECKS(0,"syntax error within surface definition"); }} CHECKS(0,"end of file encountered before end statement"); // end of file failure: // failure return *lctrptr; } /* Prints out information about all surfaces, including the surface superstructure, each surface, and panels in the surface. */ void surfaceoutput(simptr sim) { int s,p,i,dim,nident; surfacessptr srfss; surfaceptr srf; double **point,*front; printf("SURFACE PARAMETERS\n"); srfss=sim->srfss; dim=sim->dim; nident=sim->nident; if(!srfss) { printf(" No internal surfaces\n\n"); return; } printf(" Surfaces allocated: %i, surfaces defined: %i\n",srfss->maxsrf,srfss->nsrf); for(s=0;snsrf;s++) { srf=srfss->srflist[s]; printf(" Surface: %s\n",srfss->snames[s]); printf(" molecule behavior (front and back):\n"); for(i=1;iname[i],srf->faction[i],srf->baction[i]); printf(" front color: %lf %lf %lf %lf\n",srf->fcolor[0],srf->fcolor[1],srf->fcolor[2],srf->fcolor[3]); printf(" back color: %lf %lf %lf %lf\n",srf->bcolor[0],srf->bcolor[1],srf->bcolor[2],srf->bcolor[3]); printf(" edge points: %lf, polygon modes: %c %c\n",srf->edgepts,srf->fpolymode,srf->bpolymode); if(srf->maxpanel[0]) { printf(" rectangle panels allocated: %i, defined: %i\n",srf->maxpanel[0],srf->npanel[0]); for(p=0;pnpanel[0]&&p<21;p++) { if(p==20) { printf(" ...\n"); p=srf->npanel[0]-1; } point=srf->panels[0][p]->point; front=srf->panels[0][p]->front; if(dim==1) printf(" %i: %lf, facing %c0\n",p,point[0][0],front[0]==1?'+':'-'); else if(dim==2) printf(" %i: (%lf,%lf), (%lf,%lf), facing: %c%1.0f\n",p,point[0][0],point[0][1],point[1][0],point[1][1],front[0]==1?'+':'-',front[1]); else printf(" %i: (%lf,%lf,%lf), (%lf,%lf,%lf), (%lf,%lf,%lf), (%lf,%lf,%lf), facing: %c%1.0f\n",p,point[0][0],point[0][1],point[0][2],point[1][0],point[1][1],point[1][2],point[2][0],point[2][1],point[2][2],point[3][0],point[3][1],point[3][2],front[0]==1?'+':'-',front[1]); }} if(srf->maxpanel[1]) { printf(" triangle panels allocated: %i, defined: %i\n",srf->maxpanel[1],srf->npanel[1]); for(p=0;pnpanel[1]&&p<21;p++) { if(p==20) { printf(" ...\n"); p=srf->npanel[1]-1; } point=srf->panels[1][p]->point; front=srf->panels[1][p]->front; if(dim==1) printf(" %i: %lf, facing %c0\n",p,point[0][0],front[0]==1?'+':'-'); else if(dim==2) printf(" %i: (%lf,%lf), (%lf,%lf), facing: (%lf,%lf)\n",p,point[0][0],point[0][1],point[1][0],point[1][1],front[0],front[1]); else printf(" %i: (%lf,%lf,%lf), (%lf,%lf,%lf), (%lf,%lf,%lf), facing: (%lf,%lf,%lf)\n",p,point[0][0],point[0][1],point[0][2],point[1][0],point[1][1],point[1][2],point[2][0],point[2][1],point[2][2],front[0],front[1],front[2]); }} if(srf->maxpanel[2]) { printf(" sphere panels allocated: %i, defined: %i\n",srf->maxpanel[2],srf->npanel[2]); for(p=0;pnpanel[2]&&p<21;p++) { if(p==20) { printf(" ...\n"); p=srf->npanel[2]-1; } point=srf->panels[2][p]->point; front=srf->panels[2][p]->front; if(dim==1) printf(" %i: %lf, R=%lf, facing: %s\n",p,point[0][0],point[1][0],front[0]==1?"out":"in"); else if(dim==2) printf(" %i: (%lf,%lf), R=%lf, facing: %s, draw: %lf\n",p,point[0][0],point[0][1],point[1][0],front[0]==1?"out":"in",point[1][1]); else printf(" %i: (%lf,%lf,%lf), R=%lf, facing: %s, draw: %lf %lf\n",p,point[0][0],point[0][1],point[0][2],point[1][0],front[0]==1?"out":"in",point[1][1],point[1][2]); }} if(srf->maxpanel[3]) { printf(" cylinder panels allocated: %i, defined: %i\n",srf->maxpanel[3],srf->npanel[3]); for(p=0;pnpanel[3]&&p<21;p++) { if(p==20) { printf(" ...\n"); p=srf->npanel[3]-1; } point=srf->panels[3][p]->point; front=srf->panels[3][p]->front; if(dim==1) printf(" error, cylinders are not implemented in 1-D\n"); else if(dim==2) printf(" %i: (%lf,%lf) to (%lf,%lf), R=%lf, facing: %s\n",p,point[0][0],point[0][1],point[1][0],point[1][1],point[2][0],front[2]==1?"out":"in"); else printf(" %i: (%lf,%lf,%lf) to (%lf,%lf,%lf), R=%lf, facing: %s, draw: %lf %lf\n",p,point[0][0],point[0][1],point[0][2],point[1][0],point[1][1],point[1][2],point[2][0],front[0]==1?"out":"in",point[2][1],point[2][2]); }} if(srf->maxpanel[4]) { printf(" hemisphere panels allocated: %i, defined: %i\n",srf->maxpanel[4],srf->npanel[4]); for(p=0;pnpanel[4]&&p<21;p++) { if(p==20) { printf(" ...\n"); p=srf->npanel[4]-1; } point=srf->panels[4][p]->point; front=srf->panels[4][p]->front; if(dim==1) printf(" error, hemispheres are not implemented in 1-D\n"); else if(dim==2) printf(" %i: (%lf,%lf), R=%lf, facing: %s, opening: (%lf,%lf), draw: %lf\n",p,point[0][0],point[0][1],point[1][0],front[0]==1?"out":"in",point[2][0],point[2][1],point[1][1]); else printf(" %i: (%lf,%lf,%lf), R=%lf, facing: %s, opening: (%lf,%lf,%lf), draw: %lf %lf\n",p,point[0][0],point[0][1],point[0][2],point[1][0],front[0]==1?"out":"in",point[2][0],point[2][1],point[2][2],point[1][1],point[1][2]); }} printf("\n"); } return; } /* Returns the side of the panel pnl that point pt is on, which is either a ÔfÕ or a ÔbÕ for front or back, respectively. ÔbÕ is returned if the point is exactly at the panel. The value returned by this function defines the side that pt is on, so should either be called or exactly copied for other functions that care. */ char panelside(double* pt,panelptr pnl,int dim) { char face; double **point,*front,dist,cylnorm[3]; int d; point=pnl->point; front=pnl->front; if(pnl->pshape=='r') { // rectangle d=(int)front[1]; dist=front[0]*(pt[d]-point[0][d]); } else if(pnl->pshape=='t') { // triangle dist=0; for(d=0;dpshape=='s'||pnl->pshape=='h') { // sphere, hemisphere dist=0; for(d=0;dpshape=='c') { // cylinder if(dim==2) { dist=0; for(d=0;d0?'f':'b'; return face; } /* This determines if the line from pt1 to pt2 crosses the panel pnl, using a dim dimensional system. The panel includes all of its edges. If it crosses, 1 is returned, the face that is on the pt1 side of the panel is returned in faceptr, crsspt is set to the coordinates of the crossing point, and cross points to the crossing position on the line, which is a number between 0 and 1, inclusive. If it does not cross, 0 is returned and the other values are undefined. Crossing is handled very carefully such that the exact locations of pt1 and pt2, using tests that are identical to those in panelside, are used to determine which sides of the panel they are on. While crsspt will be returned with coordinates that are very close to the panel location, it may not be precisely at the panel, and there is no certainty about which side of the panel it will be on; if it matters, fix it with fixpt2panel. If the line crosses the panel more than once, which can happen for spherical or other curved panels, the smaller of the two crossing points is returned. For sphere and cylinder, 0 is returned if either both points are inside or both points are outside and the line segment does not cross the object twice. Each portion of this routine does the same things, and usually in the same order. First, the potential intersection face is determined, then the crossing value, then the crossing point, and finally it finds if intersection actually happened. For hemispheres and cylinders, if intersection does not happen for the first of two possible crossing points, it is then checked for the second point. */ int lineXpanel(double *pt1,double *pt2,panelptr pnl,char *faceptr,double *crsspt,double *cross,int dim) { surfaceptr srf; double **point,*front,dist1,dist2; double cylnorm1[3],cylnorm2[3],dot,cross2; int intsct,d; char face1,face2,faceout; srf=pnl->srf; point=pnl->point; front=pnl->front; if(pnl->pshape=='r') { // rectangle d=(int)front[1]; dist1=front[0]*(pt1[d]-point[0][d]); dist2=front[0]*(pt2[d]-point[0][d]); face1=*faceptr=dist1>0?'f':'b'; face2=dist2>0?'f':'b'; if(face1==face2) return 0; *cross=dist1/(dist1-dist2); for(d=0;dpshape=='t') { // triangle dist1=dist2=0; for(d=0;d0?'f':'b'; face2=dist2>0?'f':'b'; if(face1==face2) return 0; *cross=dist1/(dist1-dist2); for(d=0;dpshape=='s'||pnl->pshape=='h') { // sphere, hemisphere dist1=dist2=0; for(d=0;d0?'f':'b'; face2=dist2>0?'f':'b'; faceout=front[0]>0?'f':'b'; if(face1!=faceout&&face1==face2) return 0; *cross=Geo_LineXSphs(pt1,pt2,point[0],point[1][0],dim,&cross2); if(face1==face2&&!(*cross>=0&&*cross<=1)&&!(cross2>=0&&cross2<=1)) return 0; for(d=0;dpshape=='s') intsct=1; else { dot=0; for(d=0;d=0&&cross2<=1) { *cross=cross2; *faceptr=face1=='f'?'b':'f'; for(d=0;dpshape=='c') { // cylinder if(dim==2) { dist1=dist2=0; for(d=0;d0?'f':'b'; face2=dist2>0?'f':'b'; faceout=front[2]>0?'f':'b'; if(face1!=faceout&&face1==face2) return 0; *cross=Geo_LineXCyl2s(pt1,pt2,point[0],point[1],front,point[2][0],&cross2); if(face1==face2&&!(*cross>=0&&*cross<=1)&&!(cross2>=0&&cross2<=1)) return 0; } else { dist1=Geo_LineNormal3D(point[0],point[1],pt1,cylnorm1); dist2=Geo_LineNormal3D(point[0],point[1],pt2,cylnorm2); dist1=front[0]*(dist1-point[2][0]); dist2=front[0]*(dist2-point[2][0]); face1=*faceptr=dist1>0?'f':'b'; face2=dist2>0?'f':'b'; faceout=front[0]>0?'f':'b'; if(face1!=faceout&&face1==face2) return 0; *cross=Geo_LineXCyls(pt1,pt2,point[0],point[1],point[2][0],&cross2); if(face1==face2&&!(*cross>=0&&*cross<=1)&&!(cross2>=0&&cross2<=1)) return 0; } for(d=0;d=0&&cross2<=1) { *cross=cross2; *faceptr=face1=='f'?'b':'f'; for(d=0;ddim; pos1=mptr1->pos; pos2=mptr2->pos; dc1=sim->mols->difc[mptr1->ident]; dc2=sim->mols->difc[mptr2->ident]; if(dc1==0&&dc2==0) dc1=dc2=1; rxnpt=dc1/(dc1+dc2); for(bptr=mptr1->box;bptr;bptr=line2nextbox(sim,pos1,pos2,bptr)) { for(p=0;pnpanel;p++) { pnl=bptr->panel[p]; if(lineXpanel(pos1,pos2,pnl,&face,crsspt,&cross,dim)) { if(rxnptsrf->baction[mptr2->ident]:pnl->srf->faction[mptr2->ident]; else if(rxnpt>cross) action=face=='f'?pnl->srf->faction[mptr1->ident]:pnl->srf->baction[mptr1->ident]; else action=face=='f'?pnl->srf->faction[mptr1->ident]:pnl->srf->faction[mptr2->ident]; if(action!='t') return 1; }}} return 0; } /* Fixes the point pt onto the face face of panel pnl. Send in face equal to Ô0Õ if pt should be moved as close as possible to pnl. If it should also be on the front or back face of the panel, as determined by panelside, then send in face equal to ÔfÕ or ÔbÕ, respectively. This function first moves pt to the panel in a direction normal to the local panel surface and then nudges pt as required to get it to the proper side. This only considers the infinite plane of the panel, while ignoring its boundaries (similarly, hemispheres are considered to be identical to spheres and cylinders are considered to be infinitely long). */ void fixpt2panel(double *pt,panelptr pnl,int dim,char face) { int d,sign; double **point,*front,norm[3],dist,factor; point=pnl->point; front=pnl->front; dist=0; if(pnl->pshape=='r') { for(d=0;dpshape=='t') { dist=0; for(d=0;dpshape=='s'||pnl->pshape=='h') { dist=0; for(d=0;dpshape=='c') { if(dim==2) { dist=0; for(d=0;d0&&front[2]==1)||(dist<0&&front[2]==-1))?1:-1; dist=front[2]*(fabs(dist)-point[2][0]); norm[0]=sign*front[0]; norm[1]=sign*front[1]; } else if(dim==3) { dist=Geo_LineNormal3D(point[0],point[1],pt,norm); dist=front[0]*(dist-point[2][0]); }} for(d=0;dpos; front=pnl->front; if(pnl->pshape=='r') { axis=(int)front[1]; pos[axis]-=2.0*(pos[axis]-crsspt[axis]); } else if(pnl->pshape=='t') { dot=0; for(d=0;dpshape=='s'||pnl->pshape=='h') { point=pnl->point; dot=0; for(d=0;dpshape=='c') { if(dim==2) { dot=0; for(d=0;dpoint; Geo_LineNormal3D(point[0],point[1],crsspt,norm); dot=0; for(d=0;dsrf->faction[mptr->ident]; else action=pnl->srf->baction[mptr->ident]; crssptface=panelside(crsspt,pnl,dim); done=0; if(action=='t') { if(crssptface==face) fixpt2panel(crsspt,pnl,dim,face=='f'?'b':'f'); action=(face=='f')?pnl->srf->baction[mptr->ident]:pnl->srf->faction[mptr->ident]; if(action=='a') { mptr->ident=0; done=1; }} else if(action=='r') { if(crssptface!=face) fixpt2panel(crsspt,pnl,dim,face); surfacereflect(mptr,pnl,crsspt,face,dim); } else if(action=='a') { mptr->ident=0; done=1; } else if(action=='p'&&pnl->pshape=='r') { d=pnl->front[1]; srf=pnl->srf; for(p=0;pnpanel[0];p++) { pnl2=srf->panels[0][p]; if(pnl2!=pnl&&pnl2->pshape=='r'&&pnl2->front[1]==d) break; } if(pnpanel[0]) { crsspt[d]=pnl2->point[0][d]; fixpt2panel(crsspt,pnl2,dim,face); if(pnl->front[0]!=pnl2->front[0]) mptr->pos[d]=pnl2->point[0][d]+(mptr->pos[d]-pnl->point[0][d]); else mptr->pos[d]=pnl2->point[0][d]-(mptr->pos[d]-pnl->point[0][d]); if(pnl->point[0][d]point[0][d]) mptr->wrap[d]--; else mptr->wrap[d]++; }} return done; } /* Takes care of interactions between all diffusing molecules and surfaces. This transmits, reflects, or absorbs molecules, as needed, based on the panel positions and information in the molecule posx and pos elements. Absorbed molecules are killed but left in the live list with an identity of zero, for later sorting. Reflected molecules are bounced and their posx values represent the location of their last bouncing point. This function does not rely on the molecules being assigned to boxes, and nor does it assign molecules to boxes afterwards. However, it does rely on the panels being properly assigned to boxes. */ int checksurfaces(simptr sim,int ll,int reborn) { int dim,d,sctr,nmol,m,done,p,lxp,it; boxptr bptr; moleculeptr *mlist,mptr; double crossmin,crssptmin[3],crsspt[3],cross,via[3]; char face,facemin; panelptr pnl,pnlmin; if(!sim->srfss) return 0; sctr=0; dim=sim->dim; nmol=sim->mols->nl[ll]; mlist=sim->mols->live[ll]; if(reborn) m=sim->mols->topl[ll]; else m=0; for(;mposx[d]; done=0; it=0; while(!done) { if(++it>100) { for(d=0;dpos[d]=mptr->posx[d]; break; } crossmin=2; facemin='f'; pnlmin=NULL; for(bptr=pos2box(sim,via);bptr;bptr=line2nextbox(sim,via,mptr->pos,bptr)) { for(p=0;pnpanel;p++) { pnl=bptr->panel[p]; lxp=lineXpanel(via,mptr->pos,pnl,&face,crsspt,&cross,dim); if(lxp&&crossindx=NULL; bptr->nneigh=0; bptr->midneigh=0; bptr->neigh=NULL; bptr->wpneigh=NULL; bptr->nwall=0; bptr->wlist=NULL; bptr->npanel=0; bptr->panel=NULL; bptr->maxmol[0]=0; bptr->maxmol[1]=0; bptr->nmol[0]=0; bptr->nmol[1]=0; bptr->mol[0]=NULL; bptr->mol[1]=NULL; bptr->indx=allocZV(dim); if(!bptr->indx) {boxfree(bptr);return NULL;} for(d=0;dindx[d]=0; return bptr; } /* boxfree frees the box and its lists, although not the structures pointed to by the lists. */ void boxfree(boxptr bptr) { if(!bptr) return; if(bptr->indx) freeZV(bptr->indx); if(bptr->neigh) free(bptr->neigh); if(bptr->wpneigh) freeZV(bptr->wpneigh); if(bptr->wlist) free(bptr->wlist); if(bptr->panel) free(bptr->panel); if(bptr->mol[0]) free(bptr->mol[0]); if(bptr->mol[1]) free(bptr->mol[1]); free(bptr); return; } /* boxesalloc allocates and initializes an array of n boxes, including the boxes. Again, initiallization is minimal, with only the indx array of the boxes allocated, which is set to 0Õs. */ boxptr *boxesalloc(int dim,int nbox) { int b; boxptr *blist; if(nbox<1) return NULL; blist=(boxptr *) calloc(nbox,sizeof(boxptr)); if(!blist) return NULL; for(b=0;b=0;nbox--) boxfree(blist[nbox]); return; } /* boxssalloc allocates and initializes a superstructure of boxes, inculding arrays for the side, min, and size members, although the boxes are not added to the structure; i.e. blist is set to NULL and nbox is 0. Initial values for side and size members are all set to 1, min values are set to 0, and mpbox is set to 5; all of these values are typically changed later. */ boxssptr boxssalloc(int dim) { boxssptr boxs; int d; if(dim<1) return NULL; boxs=(boxssptr) malloc(sizeof(struct boxsuperstruct)); if(!boxs) return NULL; boxs->mpbox=0; boxs->boxsize=0; boxs->nbox=0; boxs->side=NULL; boxs->min=NULL; boxs->size=NULL; boxs->blist=NULL; boxs->side=allocZV(dim); if(!boxs->side) {boxssfree(boxs);return NULL;} boxs->min=(double*) calloc(dim,sizeof(double)); if(!boxs->min) {boxssfree(boxs);return NULL;} boxs->size=(double*) calloc(dim,sizeof(double)); if(!boxs->size) {boxssfree(boxs);return NULL;} for(d=0;dside[d]=1; boxs->min[d]=0; boxs->size[d]=1; } return boxs; } /* boxssfree frees a box superstructure, including the boxes. */ void boxssfree(boxssptr boxs) { if(!boxs) return; if(boxs->side) freeZV(boxs->side); if(boxs->min) free(boxs->min); if(boxs->size) free(boxs->size); if(boxs->blist) boxesfree(boxs->blist,boxs->nbox); free(boxs); return; } /* boxoutput simply lists every virtual box, along with all the details about it, where these details are the index, the number of neighbors, the neighbor mid-point, the number of maximum number of molecules, what the neighbors are, and what the wrapping codes are. As the program is currently written, this function is never called, although it could be to look for errors in box setting up. */ void boxoutput(int dim,boxssptr boxs) { int b,b2,b3,p; boxptr bptr; printf("INDIVIDUAL BOX PARAMETERS\n"); if(!boxs) { printf(" No box superstructure defined\n\n"); return; } for(b=0;bnbox;b++) { bptr=boxs->blist[b]; printf(" Box %i: indx= ",b); printZV(bptr->indx,dim); printf(" nneigh=%i midneigh=%i nwall=%i nmol,maxmol[0]=%i,%i nmol,maxmol[1]=%i,%i\n",bptr->nneigh,bptr->midneigh,bptr->nwall,bptr->nmol[0],bptr->maxmol[0],bptr->nmol[1],bptr->maxmol[1]); if(bptr->neigh) { printf(" neighbors:"); for(b2=0;b2nneigh;b2++) { b3=indx2addZV(bptr->neigh[b2]->indx,boxs->side,dim); printf(" %i",b3); } printf("\n"); } if(bptr->wpneigh) { printf(" wrap code:"); for(b2=0;b2nneigh;b2++) printf(" %i",bptr->wpneigh[b2]); printf("\n"); } printf(" %i panels",bptr->npanel); if(bptr->npanel) { printf(": "); for(p=0;pnpanel;p++) printf("%c",bptr->panel[p]->pshape); } printf("\n"); } printf("\n"); return; } /* boxssoutput displays statistics about the box superstructure, including total number of boxes, number on each side, dimensions, and the minimium position. It also prints out the requested and actual numbers of molecules per box. */ void boxssoutput(simptr sim) { int dim; boxssptr boxs; double flt1; printf("VIRTUAL BOX PARAMETERS\n"); if(!sim||!sim->boxs) { printf(" No box superstructure defined\n\n"); return; } dim=sim->dim; boxs=sim->boxs; printf(" %i boxes\n",boxs->nbox); printf(" Number of boxes on each side: "); printZV(boxs->side,dim); printf(" Minimum box position: "); printVD(boxs->min,dim); if(boxs->boxsize) printf(" Requested box width: %lf\n",boxs->boxsize); printf(" Box dimensions: "); printVD(boxs->size,dim); flt1=1.0*(sim->mols->nl[0]+sim->mols->nl[1])/boxs->nbox; if(boxs->mpbox) printf(" Requested molecules per box= %lf, actual= %lf\n",boxs->mpbox,flt1); else printf(" Molecules per box= %lf\n",flt1); printf("\n"); return; } /* expandbox is called if it turns out that a box was not allocated with enough space for molecules. bptr is a pointer to a box that needs expanding, n is the number of additional molecule spaces to add, and ll is the live list to expand. If n is negative, the box is shrunk and any molecule pointers that no longer fit are simply left out. The book keeping elements of the box are updated. The function returns 0 if it was successful and 1 if there was not enough memory for the request. */ int expandbox(boxptr bptr,int n,int ll) { moleculeptr *mlist; int m,maxmol,mn; maxmol=bptr->maxmol[ll]+n; mlist=(moleculeptr *) calloc(maxmol,sizeof(moleculeptr)); if(!mlist) return 1; mn=(n>0)?bptr->maxmol[ll]:maxmol; for(m=0;mmol[ll][m]; free(bptr->mol[ll]); bptr->mol[ll]=mlist; bptr->maxmol[ll]=maxmol; if(bptr->nmol[ll]>maxmol) bptr->nmol[ll]=maxmol; return 0;} /* pos2box returns a pointer to the box that includes the position given in pos, which is a dim size vector. If the position is outside the simulation volume, a pointer to the nearest box is returned. This routine assumes that the entire box superstructure is set up. */ boxptr pos2box(simptr sim,double *pos) { int b,d,indx,dim; boxssptr boxs; dim=sim->dim; boxs=sim->boxs; b=0; for(d=0;dmin[d])/boxs->size[d]; if(indx<0) indx=0; else if(indx>=boxs->side[d]) indx=boxs->side[d]-1; b=boxs->side[d]*b+indx; } return boxs->blist[b]; } /* Given a pointer to a box in bptr, this returns the coordinate of the low corner of the box in pos, which needs to be pre-allocated to the system dimensionality. This requires that the min and size portions of the box superstructure have been already set up. To get the other box corners, add the size values that are in the box superstructure. */ void box2pos(simptr sim,boxptr bptr,double *pos) { int d,dim; double *size,*min; dim=sim->dim; size=sim->boxs->size; min=sim->boxs->min; for(d=0;dindx[d]*size[d]; return; } /* Given a line segment which is defined by the starting point pt1 and the ending point pt2, and which is known to intersect the virtual box pointed to by bptr, this returns a pointer to the next box along the line segment. If the current box is also the final one, NULL is returned. NULL is also returned if the next box is outside of the region that is covered by virtual boxes. */ boxptr line2nextbox(simptr sim,double *pt1,double *pt2,boxptr bptr) { int dim,d,d2,side,side2,indx1,indx2,adrs; double *size,*min,crsmin,edge,crs; if(pos2box(sim,pt2)==bptr) return NULL; dim=sim->dim; size=sim->boxs->size; min=sim->boxs->min; crsmin=2; side2=0; d2=0; for(d=0;dpt1[d])?1:0; edge=min[d]+bptr->indx[d]*size[d]+(side?size[d]:0); crs=(edge-pt1[d])/(pt2[d]-pt1[d]); if(crsindx[d2]; indx2=indx1+(side2?1:-1); if(indx2<0||indx2>=sim->boxs->side[d2]) return NULL; bptr->indx[d2]=indx2; adrs=indx2addZV(bptr->indx,sim->boxs->side,dim); bptr->indx[d2]=indx1; return sim->boxs->blist[adrs]; } /* Determines if any or all of the panel pnl is in the box bptr and returns 1 if so and 0 if not. For most panel shapes, this is sufficiently complicated that this function just calls other functions in the library file Geometry.c. */ int panelinbox(simptr sim,panelptr pnl,boxptr bptr) { int dim,d,cross; double v1[3],v2[3],v3[3],v4[3],**point,*front; dim=sim->dim; box2pos(sim,bptr,v1); // v1 and v2 are set to corners of box for(d=0;dboxs->size[d]; point=pnl->point; front=pnl->front; cross=1; if(pnl->pshape=='r') { if(dim==1) cross=Geo_PtInSlab(v1,v2,point[0],dim); else if(dim==2) { v3[0]=front[1]==0?1:0; v3[1]=front[1]==0?0:1; cross=Geo_LineXaabb2(point[0],point[1],v3,v1,v2); } else if(dim==3) cross=Geo_RectXaabb3(point[0],point[1],point[3],point[0],v1,v2); } else if(pnl->pshape=='t') { if(dim==1) cross=Geo_PtInSlab(v1,v2,point[0],dim); else if(dim==2) cross=Geo_LineXaabb2(point[0],point[1],front,v1,v2); else if(dim==3) cross=Geo_TriXaabb3(point[0],point[1],point[2],front,v1,v2); } else if(pnl->pshape=='s') { if(dim==1) { if((point[0][0]-point[1][0]=v2[0])&&(point[0][0]+point[1][0]=v2[0])) cross=0; } else if(dim==2) cross=Geo_CircleXaabb2(point[0],point[1][0],v1,v2); else if(dim==3) cross=Geo_SphsXaabb3(point[0],point[1][0],v1,v2); } else if(pnl->pshape=='c') { if(dim==2) { v3[0]=point[0][0]+point[2][0]*front[0]; v3[1]=point[0][1]+point[2][0]*front[1]; v4[0]=point[1][0]+point[2][0]*front[0]; v4[1]=point[1][1]+point[2][0]*front[1]; cross=Geo_LineXaabb2(v3,v4,front,v1,v2); if(!cross) { v3[0]=point[0][0]-point[2][0]*front[0]; v3[1]=point[0][1]-point[2][0]*front[1]; v4[0]=point[1][0]-point[2][0]*front[0]; v4[1]=point[1][1]-point[2][0]*front[1]; cross=Geo_LineXaabb2(v3,v4,front,v1,v2); }} else if(dim==3) { cross=Geo_CylsXaabb3(point[0],point[1],point[2][0],v1,v2); }} else if(pnl->pshape=='h') { if(dim==2) cross=Geo_SemicXaabb2(point[0],point[1][0],point[2],v1,v2); else if(dim==3) cross=Geo_HemisXaabb3(point[0],point[1][0],point[2],v1,v2); } return cross; } /* setupboxes sets up a superstructure of boxes, and puts things in the boxes, including wall and molecule references. It requires a simulation structure with most things set up, but not box stuff; itÕs designed for the structure after itÕs returned from loadsimul. It sets up the box superstructure, then adds indicies to each box, then adds the box neighbor list along with neighbor parameters, then adds wall references to each box, and finally creates molecule lists for each box and sets both the box and molecule references to point to each other. The molecule list is 6+3Ãnmol higher than nmol to allow for more molecules; the extra pointers are all set to NULL. The function returns 0 for successful operation and 1 if it was unable to allocate sufficient memory. At the end, all simulation parameters having to do with boxes are set up. However, some lists may still be NULL, if they are empty, where these are bptr->neigh, bptr->wpneigh, and bptr->wlist. Of particular note is that bptr->wpneigh is NULL if no neighbors are wrap-around ones, for whatever reason. */ int setupboxes(simptr sim) { int dim,d,m,nbox,b,b2,w,ll; int *side,*indx; boxssptr boxs; boxptr *blist,bptr; double flt1; int *ntemp,*wptemp; int nsrf,s,ps,p,p2; surfaceptr srf; moleculeptr mptr; dim=sim->dim; boxs=sim->boxs; side=boxs->side; // box superstructure if(boxs->mpbox) { flt1=1.0; for(d=0;dwlist[2*d+1]->pos-sim->wlist[2*d]->pos; flt1=pow(1.0*(sim->mols->nd-sim->mols->topd)/boxs->mpbox/flt1,1.0/dim); } else if(boxs->boxsize) flt1=1.0/boxs->boxsize; else return 1; nbox=1; for(d=0;dmin[d]=sim->wlist[2*d]->pos; side[d]=ceil((sim->wlist[2*d+1]->pos-sim->wlist[2*d]->pos)*flt1); if(!side[d]) side[d]=1; boxs->size[d]=(sim->wlist[2*d+1]->pos-sim->wlist[2*d]->pos)/side[d]; nbox*=side[d]; } boxs->nbox=nbox; blist=boxs->blist=boxesalloc(dim,nbox); if(!blist) return 1; for(b=0;bindx,side,dim); // box->indx if(sim->accur>=3) { ntemp=allocZV(intpower(3,dim)-1); // neigh wptemp=allocZV(intpower(3,dim)-1); if(!ntemp||!wptemp) return 1; for(b=0;bindx; if(sim->accur<6) bptr->nneigh=neighborZV(indx,ntemp,side,dim,0,NULL,&bptr->midneigh); else if(sim->accur<9) { for(w=0;w<2*dim;w++) wptemp[w]=sim->wlist[w]->type=='p'; bptr->nneigh=neighborZV(indx,ntemp,side,dim,6,wptemp,&bptr->midneigh); } else { for(w=0;w<2*dim;w++) wptemp[w]=sim->wlist[w]->type=='p'; bptr->nneigh=neighborZV(indx,ntemp,side,dim,7,wptemp,&bptr->midneigh); } if(bptr->nneigh==-1) return 1; if(bptr->nneigh) { bptr->neigh=(boxptr*) calloc(bptr->nneigh,sizeof(boxptr)); if(!bptr->neigh) return 1; for(b2=0;b2nneigh;b2++) bptr->neigh[b2]=blist[ntemp[b2]]; } w=0; if(sim->accur>=6) for(b2=0;b2nneigh;b2++) w+=wptemp[b2]; if(w) { bptr->wpneigh=allocZV(bptr->nneigh); if(!bptr->wpneigh) return 1; for(b2=0;b2nneigh;b2++) bptr->wpneigh[b2]=wptemp[b2]; }} neighborZV(NULL,NULL,NULL,0,-1,NULL,NULL); freeZV(ntemp); freeZV(wptemp); } for(b=0;bnwall, wlist indx=blist[b]->indx; for(d=0;dnwall++; if(indx[d]==side[d]-1) blist[b]->nwall++; } if(blist[b]->nwall) { bptr=blist[b]; bptr->wlist=(wallptr*) calloc(bptr->nwall,sizeof(wallptr)); if(!bptr->wlist) return 1; w=0; for(d=0;dwlist[w++]=sim->wlist[2*d]; if(indx[d]==side[d]-1) bptr->wlist[w++]=sim->wlist[2*d+1]; }}} if(sim->srfss) { nsrf=sim->srfss->nsrf; for(b=0;bnpanel, panel bptr=blist[b]; for(s=0;ssrfss->srflist[s]; for(ps=0;psnpanel[ps];p++) if(panelinbox(sim,srf->panels[ps][p],bptr)) bptr->npanel++; } if(bptr->npanel) { bptr->panel=(panelptr*)calloc(bptr->npanel,sizeof(panelptr)); if(!bptr->panel) return 1; p2=0; for(s=0;ssrfss->srflist[s]; for(ps=0;psnpanel[ps];p++) if(panelinbox(sim,srf->panels[ps][p],bptr)) bptr->panel[p2++]=srf->panels[ps][p]; }}}} for(m=sim->mols->topd;mmols->nd;m++) { // box->maxmol, nmol, mol, mol->box mptr=sim->mols->dead[m]; bptr=pos2box(sim,mptr->pos); mptr->box=bptr; ll=(sim->mols->difc[mptr->ident]==0); bptr->maxmol[ll]++; } for(b=0;bmaxmol[ll]=bptr->maxmol[ll]+3*sqrt(bptr->maxmol[ll])+6; bptr->nmol[ll]=0; bptr->mol[ll]=(moleculeptr*) calloc(bptr->maxmol[ll],sizeof(moleculeptr)); if(!bptr->mol[ll]) return 1; for(m=0;mmaxmol[ll];m++) bptr->mol[ll][m]=NULL; }} return 0; } /* assignmolecs puts molecules in boxes by overwriting the lists of molecules that are in each box with molecules from mols. It only assigns the live list number ll. Molecules that are outside the set of boxes are assigned to the nearest box. If more molecules belong in a box than actually fit, 5 more spaces are allocated using expandbox. The function returns 0 unless memory could not be allocated by expandbox, in which case only some of the molecules are assigned and it returns 1. */ int assignmolecs(simptr sim,int ll) { int m,nmol,b; boxptr *blist,bptr; boxssptr boxs; moleculeptr *mlist; boxs=sim->boxs; blist=boxs->blist; mlist=sim->mols->live[ll]; nmol=sim->mols->nl[ll]; for(b=0;bnbox;b++) blist[b]->nmol[ll]=0; for(m=0;mpos); if(bptr->nmol[ll]==bptr->maxmol[ll]) if(expandbox(bptr,5,ll)) return 1; bptr->mol[ll][bptr->nmol[ll]++]=mlist[m]; mlist[m]->box=bptr; } return 0; } /* Updates the information about which molecule is in which box, for those molecules that are in master list ll (0 for diffusing, 1 for not). If reborn is 0, all molecules in list ll are updated; otherwise, only the reborn ones are. This assumes that all molecules that are being checked were in a box, meaning that the box element of the molecule structure listed a box and that the mol list of that box listed the molecule. Molecules are arranged in boxes according to the location of the pos element of the molecules. Molecules outside the set of boxes are assigned to the nearest box. If more molecules belong in a box than actually fit, the number of spaces is doubled using expandbox. The function returns 0 unless memory could not be allocated by expandbox, in which case it fails and returns 1. */ int reassignmolecs(simptr sim,int ll,int reborn) { int m,nmol,m2; boxptr bptr; moleculeptr mptr,*mlist,*mlist2; mlist=sim->mols->live[ll]; nmol=sim->mols->nl[ll]; if(reborn) { mlist+=sim->mols->topl[ll]; nmol-=sim->mols->topl[ll]; } for(m=0;mpos); if(mptr->box!=bptr) { mlist2=mptr->box->mol[ll]; for(m2=0;mlist2[m2]!=mptr;m2++); mlist2[m2]=mlist2[--mptr->box->nmol[ll]]; mptr->box=bptr; if(bptr->nmol[ll]==bptr->maxmol[ll]) if(expandbox(bptr,bptr->nmol[ll],ll)) return 1; bptr->mol[ll][bptr->nmol[ll]++]=mptr; }} return 0; } /******************************************************************************/ /***************************** Simulation structure ***************************/ /******************************************************************************/ /* simalloc allocates a simulation structure. The difc and difm lists are allocated and initialized. Default diffusion matrics are all 0Õs, except with Ð1 in the first element. Walls are allocated and inialized. The box superstructure is allocated and initialized, although the list of boxes is left as NULL. The molecule superstructure is left as NULL. The commands superstructure is allocated, but the queue of commands and the output file lists are left as NULLs. root is required for the command superstructure. */ simptr simalloc(int dim,int maxident,char *root) { simptr sim; int i,max; sim=NULL; CHECK(dim>0&&maxident>0); CHECK(sim=(simptr) malloc(sizeof(struct simstruct))); sim->randseed=randomize(); sim->dim=dim; sim->maxident=maxident; sim->nident=1; sim->name=NULL; sim->graphics=0; sim->graphicit=1; sim->tiffit=0; sim->framepts=2; sim->gridpts=0; sim->framecolor[0]=0; sim->framecolor[1]=0; sim->framecolor[2]=0; sim->framecolor[3]=1; sim->backcolor[0]=1; sim->backcolor[1]=1; sim->backcolor[2]=1; sim->backcolor[3]=1; sim->accur=10; sim->time=0; sim->tmin=0; sim->tmax=10; sim->dt=1; sim->rxn[0]=sim->rxn[1]=sim->rxn[2]=NULL; sim->mols=NULL; sim->wlist=NULL; sim->srfss=NULL; sim->boxs=NULL; sim->cmds=NULL; CHECK(sim->name=(char **) calloc(maxident,sizeof(char *))); for(i=0;iname[i]=NULL; for(i=0;iname[i]=EmptyString()); } strncpy(sim->name[0],"empty",STRCHAR); CHECK(sim->wlist=wallsalloc(dim)); CHECK(sim->boxs=boxssalloc(dim)); CHECK(sim->cmds=scmdssalloc(&docommand,(void*)sim,root)); max=dim>maxident?dim:maxident; CHECK(sim->v1=(double*) calloc(max,sizeof(double))); CHECK(sim->v2=(double*) calloc(max,sizeof(double))); CHECK(sim->v3=(double*) calloc(max,sizeof(double))); CHECK(sim->m1=(double*) calloc(dim*dim,sizeof(double))); CHECK(sim->m2=(double*) calloc(dim*dim,sizeof(double))); CHECK(sim->m3=(double*) calloc(dim*dim,sizeof(double))); CHECK(sim->z1=(int*) calloc(max,sizeof(int))); CHECK(sim->z2=(int*) calloc(max,sizeof(int))); CHECK(sim->z3=(int*) calloc(max,sizeof(int))); for(i=0;iv1[i]=sim->v2[i]=sim->v3[i]=0; for(i=0;im1[i]=sim->m2[i]=sim->m3[i]=0; for(i=0;iz1[i]=sim->z2[i]=sim->z3[i]=0; return sim; failure: simfree(sim); return NULL; } /* simfree frees a simulation strucutre, including every part of everything in it. */ void simfree(simptr sim) { int i,dim; int maxident; if(!sim) return; dim=sim->dim; maxident=sim->maxident; if(sim->rxn[0]) rxnfree(sim->rxn[0],maxident); if(sim->rxn[1]) rxnfree(sim->rxn[1],maxident); if(sim->rxn[2]) rxnfree(sim->rxn[2],maxident); sim->rxn[0]=sim->rxn[1]=sim->rxn[2]=NULL; if(sim->mols) molssfree(sim->mols,maxident); sim->mols=NULL; if(sim->wlist) wallsfree(sim->wlist,dim); sim->wlist=NULL; if(sim->srfss) surfacessfree(sim->srfss); sim->srfss=NULL; if(sim->boxs) boxssfree(sim->boxs); sim->boxs=NULL; if(sim->cmds) scmdssfree(sim->cmds); sim->cmds=NULL; if(maxident&&sim->name) for(i=0;iname[i]) free(sim->name[i]); if(sim->name) free(sim->name); if(sim->v1) free(sim->v1); if(sim->v2) free(sim->v2); if(sim->v3) free(sim->v3); if(sim->m1) freeM(sim->m1); if(sim->m2) freeM(sim->m2); if(sim->m3) freeM(sim->m3); if(sim->z1) free(sim->z1); if(sim->z2) free(sim->z2); if(sim->z3) free(sim->z3); free(sim); return; } /* loadsimul loads all simulation parameters from a configuration file, using a format described above. fileroot is sent in as the root of the filename, including all colons, slashes, or backslashes; if the configuration file is in the same directory as Smoldyn, fileroot should be an empty string. filename is sent in as just the file name and any extension. erstr is sent in as an empty string of size STRCHAR and is returned with an error message if an error occurs. smptr is sent in as a pointer to the variable that will point to the simstruct; it is returned pointing to a pointer to an initiallized simstruct. This routine calls loadrxn to load in any reactions. The following things are set up after this routine is completed: all molecule elements except box; all molecule superstructure elements; all wall elements; box superstructure element mpbox, but no other elements; no boxes are allocated or set up; all reaction structure elements except rate2 and the product template position vectors (pos in each product); the command superstructure, including all of its elements; and all simulation structure elements except for sub-elements that have already been listed. If the configuration file loads successfully, the routine returns 0. If the file could not be found, it returns 10 and an error message. If an error was caught during file loading, the return value is 10 plus the line number of the file with an error, along with an error message. If there is an error, all structures are freed automatically. */ int loadsimul(simptr *smptr,char *fileroot,char *filename,char *erstr) { int filenum; struct {int lctr;char fname[STRCHAR];FILE *fptr; } filestack[32]; int done,itct,got[18],er; char word[STRCHAR],line[STRCHAR],nm[STRCHAR],*line2,*linetest,ch; char **name; simptr sim; molssptr mols; moleculeptr *dead; wallptr *wlist; cmdssptr cmds; int dim,maxident,nident,graphics,graphicit,tiffit,incomment; int i,m,nmol,d,topd; int i1; double flt1,flt2,accur; topd=0; dead=NULL; mols=NULL; name=NULL; cmds=NULL; wlist=NULL; nident=0; if(!erstr||!filename||!fileroot) return 0; *smptr=sim=NULL; filenum=0; done=0; incomment=0; setstdZV(got,18,0); graphics=0; graphicit=1; tiffit=0; accur=10; filestack[0].lctr=10; strncpy(filestack[0].fname,fileroot,STRCHAR); strncat(filestack[0].fname,filename,STRCHAR-strlen(fileroot)); filestack[0].fptr=fopen(filestack[0].fname,"r"); CHECKS(filestack[0].fptr,"file not found"); while(!done) { linetest=fgets(line,STRCHAR,filestack[filenum].fptr); if(linetest) { filestack[filenum].lctr++; if(strchr(line,'#')) *strchr(line,'#')='\0';// comment itct=sscanf(line,"%s",word); line2=strnword(line,2); if(line2&&strchr(line2,'\n')) *(strchr(line2,'\n'))='\0'; } if(!linetest||(itct&&!incomment&&!strcmp(word,"end_file"))) {// end_file fclose(filestack[filenum].fptr); if(--filenum<0) done=1; } else if(itct<=0); // blank line else if(!strcmp(word,"/*")) { // start comment with /* incomment=1; } else if(!strcmp(word,"*/")) { // end comment with */ CHECKS(incomment,"*/ without a corresponding /*"); incomment=0; } else if(incomment); else if(!strcmp(word,"start_reaction")) { // start_reaction CHECKS(got[1],"need to enter names before reactions"); CHECKS(!line2,"unexpected text following start_reaction"); CHECKS(!loadrxn(sim,filestack[filenum].fptr,&(filestack[filenum].lctr),erstr),erstr); } else if(!strcmp(word,"start_surface")) { // start_surface CHECKS(got[17],"need to enter max_surface before start_surface"); CHECKS(!line2,"unexpected text following start_surface"); CHECKS(!loadsurface(sim,filestack[filenum].fptr,&(filestack[filenum].lctr),erstr),erstr); } else if(!line2) { // just word CHECKS(0,"unknown word or missing parameter"); } else if(!strcmp(word,"read_file")) { // read_file CHECKS(filenum<31&&filenum0,"dim needs to be >0"); CHECKS(!strnword(line2,2),"unexpected text following dim"); } else if(!strcmp(word,"max_names")) { // max_names, got[1] CHECKS(got[0],"dim has to be entered before max_names"); CHECKS(!got[1],"max_names can only be entered once, and not after 'names'"); got[1]++; itct=sscanf(line2,"%i",&maxident); CHECKS(itct==1,"max_names needs to be an integer"); CHECKS(maxident>0,"max_names needs to be >0"); nident=1; CHECKS(sim=*smptr=simalloc(dim,maxident,fileroot),"out of memory in simalloc()"); sim->nident=nident; wlist=sim->wlist; cmds=sim->cmds; name=sim->name; CHECKS(!strnword(line2,2),"unexpected text following max_names"); } else if(!strcmp(word,"name")) { // name CHECKS(got[0],"dim has to be entered before name"); CHECKS(got[1],"max_names has to be entered before name"); CHECKS(nidentnident=nident; CHECKS(itct==1,"error reading molecule name"); CHECKS(!strnword(line2,2),"unexpected text following name"); } else if(!strcmp(word,"names")) { // names, got[1] CHECKS(got[0],"dim has to be entered before names"); CHECKS(!got[1],"names can only be entered once, and not after 'max_names'"); got[1]++; nident=1+wordcount(line2); maxident=nident; CHECKS(!strstr(line2,"empty"),"'empty' is not a permitted name"); CHECKS(sim=*smptr=simalloc(dim,maxident,fileroot),"out of memory in simalloc()"); sim->nident=nident; wlist=sim->wlist; cmds=sim->cmds; name=sim->name; for(i=1;i=0,"max_mol needs to be >0"); CHECKS(sim->mols=mols=molssalloc(dim,i1,maxident),"out of memory in molssalloc()"); dead=mols->dead; topd=mols->topd; CHECKS(!strnword(line2,2),"unexpected text following max_mol"); } else if(!strcmp(word,"molperbox")) { // molperbox, got[5] CHECKS(got[1],"need to enter names before molperbox"); CHECKS(!got[5],"molperbox or boxsize can only be entered once"); got[5]++; itct=sscanf(line2,"%lf",&sim->boxs->mpbox); CHECKS(itct==1,"molperbox needs to be a number"); CHECKS(sim->boxs->mpbox>0,"molperbox needs to be >0"); CHECKS(!strnword(line2,2),"unexpected text following molperbox"); } else if(!strcmp(word,"boxsize")) { // boxsize, got[5] CHECKS(got[1],"need to enter names before boxsize"); CHECKS(!got[5],"molperbox or boxsize can only be entered once"); got[5]++; itct=sscanf(line2,"%lf",&sim->boxs->boxsize); CHECKS(itct==1,"boxsize needs to be a number"); CHECKS(sim->boxs->boxsize>0,"boxsize needs to be >0"); CHECKS(!strnword(line2,2),"unexpected text following boxsize"); } else if(!strcmp(word,"difc")) { // difc CHECKS(got[3],"need to enter max_mol before difc"); itct=sscanf(line2,"%s %lf",nm,&flt1); CHECKS(itct==2,"difc format: name value"); i=stringfind(name,nident,nm); CHECKS(i>0,"name not recognized for difc"); CHECKS(flt1>=0,"diffusion constant needs to be >=0"); mols->difc[i]=flt1; CHECKS(!strnword(line2,3),"unexpected text following difc"); } else if(!strcmp(word,"difm")) { // difm CHECKS(got[3],"need to enter max_mol before difm"); itct=sscanf(line2,"%s",nm); CHECKS(itct==1,"difm format: name matrix"); i=stringfind(name,nident,nm); CHECKS(i>0,"name not recognized for difm") CHECKS(line2=strnword(line2,2),"missing matrix in difm"); if(!mols->difm[i]) { CHECKS(mols->difm[i]=(double*)calloc(dim*dim,sizeof(double)),"out of memory in difm"); } itct=strreadnd(line2,dim*dim,mols->difm[i],NULL); CHECKS(itct==dim*dim,"incomplete matrix in difm"); CHECKS(!strnword(line2,dim*dim+1),"unexpected text following difm"); } else if(!strcmp(word,"color")) { // color CHECKS(got[3],"need to enter max_mol before color"); itct=sscanf(line2,"%s",nm); CHECKS(itct==1,"color format: name red green blue"); i=stringfind(name,nident,nm); CHECKS(i>0,"name not recognized for color"); CHECKS(line2=strnword(line2,2),"missing color data for color"); itct=sscanf(line2,"%lf %lf %lf",&(mols->color[i][0]),&(mols->color[i][1]),&(mols->color[i][2])); CHECKS(itct==3,"color format: name red green blue"); CHECKS(mols->color[i][0]>=0&&mols->color[i][0]<=1,"color value needs to between 0 and 1"); CHECKS(mols->color[i][1]>=0&&mols->color[i][1]<=1,"color value needs to between 0 and 1"); CHECKS(mols->color[i][2]>=0&&mols->color[i][2]<=1,"color value needs to between 0 and 1"); CHECKS(!strnword(line2,4),"unexpected text following color"); } else if(!strcmp(word,"display_size")) { // display_size CHECKS(got[3],"need to enter max_mol before display_size"); itct=sscanf(line2,"%s",nm); CHECKS(itct==1,"display_size format: name size"); i=stringfind(name,nident,nm); CHECKS(i>0,"name not recognized for display_size"); CHECKS(line2=strnword(line2,2),"missing data for display_size"); itct=sscanf(line2,"%lf",&(mols->display[i])); CHECKS(itct==1,"display_size format: name size"); CHECKS(!strnword(line2,2),"unexpected text following display_size"); } else if(!strcmp(word,"time_start")) { // time_start, got[6] CHECKS(got[1],"need to enter names before time_start"); CHECKS(!got[6],"time_start can only be entered once"); got[6]++; itct=sscanf(line2,"%lf",&sim->tmin); CHECKS(itct==1,"time_start needs to be a number"); CHECKS(!strnword(line2,2),"unexpected text following time_start"); } else if(!strcmp(word,"time_stop")) { // time_stop, got[7] CHECKS(got[1],"need to enter names before time_stop"); CHECKS(!got[7],"time_stop can only be entered once"); got[7]++; itct=sscanf(line2,"%lf",&sim->tmax); CHECKS(itct==1,"time_stop needs to be a number"); CHECKS(!strnword(line2,2),"unexpected text following time_stop"); } else if(!strcmp(word,"time_step")) { // time_step, got[8] CHECKS(got[1],"need to enter names before time_step"); CHECKS(!got[8],"time_step can only be entered once"); got[8]++; itct=sscanf(line2,"%lf",&sim->dt); CHECKS(itct==1,"time_step needs to be a number"); CHECKS(sim->dt>0,"time step must be >0"); CHECKS(!strnword(line2,2),"unexpected text following time_step"); } else if(!strcmp(word,"time_now")) { // time_now, got[2] CHECKS(got[1],"need to enter names before time_now"); CHECKS(!got[2],"time_now can only be entered once"); got[2]++; itct=sscanf(line2,"%lf",&sim->time); CHECKS(itct==1,"time_now needs to be a number"); CHECKS(!strnword(line2,2),"unexpected text following time_now"); } else if(!strcmp(word,"rand_seed")) { // rand_seed, got[9] CHECKS(got[1],"need to enter names before rand_seed"); CHECKS(!got[9],"rand_seed can only be entered once"); got[9]++; itct=sscanf(line2,"%i",&i1); CHECKS(itct==1,"rand_seed needs to be an integer"); sim->randseed=i1; srand(i1); CHECKS(!strnword(line2,2),"unexpected text following rand_seed"); } else if(!strcmp(word,"graphics")) { // graphics, got[10] CHECKS(!got[10],"graphics can only be entered once"); got[10]++; itct=sscanf(line2,"%s",nm); CHECKS(itct==1,"misformated graphics command"); if(!strcmp(nm,"none")) graphics=0; else if(!strcmp(nm,"opengl")) graphics=1; else if(!strcmp(nm,"opengl_good")) graphics=2; else CHECKS(0,"graphics method not recognized"); CHECKS(!strnword(line2,2),"unexpected text following graphics"); } else if(!strcmp(word,"graphic_iter")) { // graphic_iter, got[11] CHECKS(!got[11],"graphic_iter can only be entered once"); got[11]++; itct=sscanf(line2,"%i",&graphicit); CHECKS(itct==1,"graphic_iter need to be a number"); CHECKS(graphicit>0,"graphicit needs to be >=1"); CHECKS(!strnword(line2,2),"unexpected text following graphic_iter"); } else if(!strcmp(word,"tiff_iter")) { // tiff_iter, got[12] CHECKS(!got[12],"tiff_iter can only be entered once"); got[12]++; itct=sscanf(line2,"%i",&tiffit); CHECKS(itct==1,"tiff_iter needs to be a number"); CHECKS(!strnword(line2,2),"unexpected text following tiff_iter"); } else if(!strcmp(word,"tiff_name")) { // tiff_name, got[16] CHECKS(!got[16],"tiff_name can only be entered once"); got[16]++; itct=sscanf(line2,"%s",nm); CHECKS(itct==1,"tiff_name needs to be a string"); gl2SetOptionStr("TiffName",nm); CHECKS(!strnword(line2,2),"unexpected text following tiff_name"); } else if(!strcmp(word,"tiff_min")) { // tiff_min, got[14] CHECKS(!got[14],"tiff_min can only be entered once"); got[14]++; itct=sscanf(line2,"%i",&i1); CHECKS(itct==1,"tiff_min needs to be a number"); gl2SetOptionInt("TiffNumber",i1); CHECKS(!strnword(line2,2),"unexpected text following tiff_min"); } else if(!strcmp(word,"tiff_max")) { // tiff_max, got[15] CHECKS(!got[15],"tiff_max can only be entered once"); got[15]++; itct=sscanf(line2,"%i",&i1); CHECKS(itct==1,"tiff_max needs to be a number"); gl2SetOptionInt("TiffNumMax",i1); CHECKS(!strnword(line2,2),"unexpected text following tiff_max"); } else if(!strcmp(word,"frame_thickness")) { // frame_thickness CHECKS(got[1],"need to enter names before frame_thickness"); itct=sscanf(line2,"%lf",&sim->framepts); CHECKS(itct==1,"frame_thickness needs to be a number"); CHECKS(sim->framepts>=0,"frame_thickness needs to be ³0"); CHECKS(!strnword(line2,2),"unexpected text following frame_thickness"); } else if(!strcmp(word,"grid_thickness")) { // grid_thickness CHECKS(got[1],"need to enter names before grid_thickness"); itct=sscanf(line2,"%lf",&sim->gridpts); CHECKS(itct==1,"grid_thickness needs to be a number"); CHECKS(sim->gridpts>=0,"grid_thickness needs to be ³0"); CHECKS(!strnword(line2,2),"unexpected text following grid_thickness"); } else if(!strcmp(word,"frame_color")) { // frame_color CHECKS(got[1],"need to enter names before frame_color"); itct=sscanf(line2,"%lf %lf %lf",&(sim->framecolor[0]),&(sim->framecolor[1]),&(sim->framecolor[2])); CHECKS(itct==3,"frame_color format: red green blue [alpha]"); CHECKS(sim->framecolor[0]>=0&&sim->framecolor[0]<=1,"frame_color values need to between 0 and 1"); CHECKS(sim->framecolor[1]>=0&&sim->framecolor[1]<=1,"frame_color values need to between 0 and 1"); CHECKS(sim->framecolor[2]>=0&&sim->framecolor[2]<=1,"frame_color values need to between 0 and 1"); if((line2=strnword(line2,4))!=NULL) { itct=sscanf(line2,"%lf",&(sim->framecolor[3])); CHECKS(sim->framecolor[3]>=0&&sim->framecolor[3]<=1,"frame_color alpha value needs to be between 0 and 1"); line2=strnword(line2,2); } CHECKS(!line2,"unexpected text following frame_color"); } else if(!strcmp(word,"background_color")) { // background_color CHECKS(got[1],"need to enter names before background_color"); itct=sscanf(line2,"%lf %lf %lf",&(sim->backcolor[0]),&(sim->backcolor[1]),&(sim->backcolor[2])); CHECKS(itct==3,"background_color format: red green blue [alpha]"); CHECKS(sim->backcolor[0]>=0&&sim->backcolor[0]<=1,"background_color values need to between 0 and 1"); CHECKS(sim->backcolor[1]>=0&&sim->backcolor[1]<=1,"background_color values need to between 0 and 1"); CHECKS(sim->backcolor[2]>=0&&sim->backcolor[2]<=1,"background_color values need to between 0 and 1"); if((line2=strnword(line2,4))!=NULL) { itct=sscanf(line2,"%lf",&(sim->backcolor[3])); CHECKS(sim->backcolor[3]>=0&&sim->backcolor[3]<=1,"background_color alpha value needs to be between 0 and 1"); line2=strnword(line2,2); } CHECKS(!line2,"unexpected text following background_color"); } else if(!strcmp(word,"accuracy")) { // accuracy, got[13] CHECKS(!got[13],"accuracy can only be entered once"); got[13]++; itct=sscanf(line2,"%lf",&accur); CHECKS(itct==1,"accuracy needs to be a number"); CHECKS(!strnword(line2,2),"unexpected text following accuracy"); } else if(!strcmp(word,"low_wall")) { // low_wall CHECKS(got[1],"need to enter names before low_wall"); itct=sscanf(line2,"%i %lf %c",&d,&flt1,&ch); CHECKS(itct==3,"low_wall format: dimension position type"); CHECKS(d>=0&&dpos=flt1; CHECKS(ch=='r'||ch=='p'||ch=='a'||ch=='t',"low_wall type needs to be r, p, a, or t"); wlist[2*d]->type=ch; CHECKS(!strnword(line2,4),"unexpected text following low_wall"); } else if(!strcmp(word,"high_wall")) { // high_wall CHECKS(got[1],"need to enter names before high_wall"); itct=sscanf(line2,"%i %lf %c",&d,&flt1,&ch); CHECKS(itct==3,"high_wall format: dimension position type"); CHECKS(d>=0&&dpos=flt1; CHECKS(ch=='r'||ch=='p'||ch=='a'||ch=='t',"high_wall type needs to be r, p, a, or t"); wlist[2*d+1]->type=ch; CHECKS(!strnword(line2,4),"unexpected text following high_wall"); } else if(!strcmp(word,"max_surface")) { // max_surface, got[17] CHECKS(got[1],"need to enter names before max_surface"); CHECKS(!got[17],"max_surface can only be entered once"); got[17]++; itct=sscanf(line2,"%i",&i1); CHECKS(itct==1,"max_surface needs to be a number"); CHECKS(i1>=0,"max_surface must be at least 0"); if(i1) { sim->srfss=surfacessalloc(i1,maxident,dim); CHECKS(sim->srfss,"memory allocation failure at max_surface"); } CHECKS(!strnword(line2,2),"unexpected text following max_surface"); } else if(!strcmp(word,"mol")) { // mol CHECKS(got[3],"need to enter max_mol before mol"); itct=sscanf(line2,"%i %s",&nmol,nm); CHECKS(itct==2,"mol format: number name position_vector"); CHECKS(nmol>=0,"number of molecules added needs to be >=0"); CHECKS(topd-nmol>=0,"more molecules assigned than allocated"); i=stringfind(name,nident,nm); CHECKS(i>0,"name not recognized for mol"); for(m=topd-1;m>=topd-nmol;m--) dead[m]->ident=i; CHECKS(line2=strnword(line2,2),"insufficient data in mol command"); for(d=0;d=topd-nmol;m--) dead[m]->pos[d]=dead[m]->posx[d]=unirand30(wlist[2*d]->pos,wlist[2*d+1]->pos); } else if(strchr(line2+1,'-')&&(!strchr(line2,' ')||strchr(line2+1,'-')=topd-nmol;m--) dead[m]->pos[d]=dead[m]->posx[d]=unirand30(flt1,flt2); } else { itct=sscanf(line2,"%lf",&flt1); CHECKS(itct==1,"cannot read position value for mol"); for(m=topd-1;m>=topd-nmol;m--) dead[m]->pos[d]=dead[m]->posx[d]=flt1; }} topd-=nmol; mols->topd=topd; CHECKS(!strnword(line2,2),"unexpected text following mol"); } else if(!strcmp(word,"output_root")) { // output_root CHECKS(got[1],"names have to be entered before output_root"); CHECKS(!scmdsetfroot(cmds,line2),"output_root can only be entered once"); } else if(!strcmp(word,"output_files")) { // output_files CHECKS(got[1],"need to enter names before output_files"); er=scmdsetfnames(cmds,line2); CHECKS(er!=1,"out of memory in output_files"); CHECKS(er!=2,"error reading file name"); CHECKS(er!=3,"output_files can only be entered once"); CHECKS(er!=4,"BUG: variable cmds became NULL"); } else if(!strcmp(word,"output_file_number")) { // output_file_number CHECKS(got[1],"need to enter names before output_file_number"); itct=sscanf(line2,"%s %i",nm,&i1); CHECKS(itct==2,"format for output_file_number: filename number"); CHECKS(i1>=0,"output_file_number needs to be ³0"); CHECKS(!scmdsetfsuffix(cmds,nm,i1),"error setting output_file_number"); CHECKS(!strnword(line2,3),"unexpected text following output_file_number"); } else if(!strcmp(word,"max_cmd")) { // max_cmd, got[4] CHECKS(got[1],"need to enter names before max_cmd"); CHECKS(!got[4],"max_cmd can only be entered once and has to be before any commands"); got[4]++; itct=sscanf(line2,"%i",&i1); CHECKS(itct==1,"max_cmd needs to be an integer"); CHECKS(i1>=0,"max_cmd needs to be >=0"); CHECKS(!scmdqalloc(cmds,i1),"command queue allocation error"); CHECKS(!strnword(line2,2),"unexpected text following max_cmd"); } else if(!strcmp(word,"cmd")) { // cmd CHECKS(got[1],"need to enter names before commands"); CHECKS(got[6],"need to enter time_start before commands"); CHECKS(got[7],"need to enter time_stop before commands"); CHECKS(got[8],"need to enter time_step before commands"); if(!got[4]) got[4]++; er=scmdstr2cmd(cmds,line2,sim->tmin,sim->tmax,sim->dt); CHECKS(er!=1,"out of memory in cmd"); CHECKS(er!=2,"BUG: no command superstructure for cmd"); CHECKS(er!=3,"cmd format: type [on off dt] string"); CHECKS(er!=4,"command string is missing"); CHECKS(er!=5,"cmd time step needs to be >0"); CHECKS(er!=6,"command timing type character not recognized"); CHECKS(er!=7,"insertion of command in queue failed"); CHECKS(er!=8,"cmd time multiplier needs to be >1"); } else { // unknown word CHECKS(0,"syntax error"); }} sim->graphics=graphics; sim->graphicit=graphicit; sim->tiffit=tiffit; sim->accur=accur; CHECKS(got[0],"dimensionality was not found in configuration file"); // assorted error checking CHECKS(got[1],"names were not found in configuration file"); CHECKS(got[3],"no molecules were listed in configuration file"); CHECKS(got[6],"no start time was listed in configuration file"); CHECKS(got[7],"no stop time was listed in configuration file"); CHECKS(got[8],"no time step was listed in configuration file"); if(!got[2]) sim->time=sim->tmin; if(!got[5]) sim->boxs->mpbox=5; for(d=0;dpospos,"low_wall positions need to be smaller than high_wall positions"); *smptr=sim; return 0; failure: // failure if(done) i1=filestack[0].lctr+1; else { i1=filestack[filenum].lctr; strncat(erstr,", file ",STRCHAR-strlen(erstr)); strncat(erstr,filestack[filenum].fname,STRCHAR-strlen(erstr)); for(;filenum>=0;filenum--) if(filestack[filenum].fptr) fclose(filestack[filenum].fptr); } simfree(sim); return i1; } /* setupstructs sets up and loads values for all the structures as well as global variables. This routine calls the other initialization routines, so they do not have to be called from elsewhere. Other minor things are set up here, including setting the lookup table for normally distributed random numbers. It also displays the status to stdout and calls output routines for each structure, allowing verification of the initiallization. Send in root and name with strings for the path and name of the input file. vb is a flag for verbose operation, 1 for verbose and 0 for quiet. It returns 0 for correct operation and 1 for an error. If it succeeds, smptr is returned pointing to a simulation structure. Otherwise, smptr is set to NULL and an error messages is displayed on stderr. */ int setupstructs(char *root,char *name,simptr *smptr,int vb) { simptr sim; int er,order,dim,r; char erstr[STRCHAR]; sim=*smptr=NULL; if(vb&&name) { printf("\nCONFIGURATION FILE\n"); printf(" Root: '%s'\n",root); printf(" Name: '%s'\n",name); } er=loadsimul(&sim,root,name,erstr); if(er) { fprintf(stderr," Error reading file in line %i\n",er-10); fprintf(stderr," %s\n",erstr); return 1; } if(vb) printf(" Loaded file successfully\n\n"); dim=sim->dim; randtableD(GaussTable,RANDTABLEMAX+1,1); randshuffletableD(GaussTable,RANDTABLEMAX+1); setdiffusion(sim); er=setupboxes(sim); if(er) { fprintf(stderr,"out of memory setting up spatial partitions\n"); simfree(sim); return 1; } er=molsort(sim->mols,0); if(er) { fprintf(stderr,"out of memory during molecule sorting\n"); simfree(sim); return 1; } if(vb) simoutput(sim); if(vb) walloutput(dim,sim->wlist); if(vb) molssoutput(sim); if(vb) surfaceoutput(sim); if(vb) scmdoutput(sim->cmds); if(vb) boxssoutput(sim); // if(vb) boxoutput(dim,sim->boxs); // This line only used for debugging for(order=0;order<3;order++) { for(r=0;sim->rxn[order]&&rrxn[order]->total;r++) if(sim->rxn[order]->rate[r]>=0&&sim->rxn[order]->rate2[r]>=0&&vb) printf("WARNING: rate and rate_internal are defined for reaction order %i, name %s\n",order,sim->rxn[order]->rname[r]); er=setrates(sim,order); if(er>=0) { fprintf(stderr,"Error setting rate for reaction order %i, reaction %s\n",order,sim->rxn[order]->rname[er]); fprintf(stderr," (check product parameter if it is reversible, and that a reactant diffuses.)\n"); simfree(sim); return 1; }} for(order=0;order<3;order++) { erstr[0]='\0'; er=setproducts(sim,order,erstr); if(er>=0) { fprintf(stderr,"Error setting products for reaction order %i, reaction %s\n",order,sim->rxn[order]->rname[er]); fprintf(stderr,"%s\n",erstr); simfree(sim); return 1; } if(strlen(erstr)) fprintf(stderr,"%s\n",erstr); } for(order=0;order<3;order++) if(vb&&sim->rxn[order]) rxnoutput(sim,order); if(vb) checkparams(sim); *smptr=sim; return 0; } /* simoutput prints out the overall simulation parameters, including simulation time information, graphics information, the number of dimensions, what the molecule types are, the output files, and the accuracy. */ void simoutput(simptr sim) { int i,dim,nident; printf("SIMULATION PARAMETERS\n"); if(!sim) { printf(" No simulation parameters\n\n"); return; } dim=sim->dim; nident=sim->nident; printf(" %i dimensions\n",dim); printf(" Accuracy level: %lf\n",sim->accur); printf(" Random number seed: %u\n",sim->randseed); printf(" %i identities (excluding empty molecules) of %i allocated:\n",sim->nident-1,sim->maxident); printf(" %s",sim->name[1]); for(i=2;iname[i]); printf("\n"); printf(" Time from %lf to %lf step %lf\n",sim->tmin,sim->tmax,sim->dt); if(sim->time!=sim->tmin) printf(" Current time: %lf\n",sim->time); printf(" Graphical output: %s",sim->graphics?"OpenGL":"none"); if(sim->graphics) printf(", every %i iterations\n",sim->graphicit); else printf("\n"); if(sim->graphics) { printf(" frame thickness: %lf",sim->framepts); if(sim->gridpts) printf(", grid thickness: %lf",sim->gridpts); printf("\n"); printf(" frame color: %lf,%lf,%lf,%lf\n",sim->framecolor[0],sim->framecolor[1],sim->framecolor[2],sim->framecolor[3]); printf(" background color: %lf,%lf,%lf,%lf\n",sim->backcolor[0],sim->backcolor[1],sim->backcolor[2],sim->backcolor[3]); } printf("\n"); return; } /* checkparams checks that the simulation parameters, including parameters of sub-structures, have reasonable values. If values seem to be too small or too large, a warning is displayed to the standard output, although this does not affect continuation of the program. */ void checkparams(simptr sim) { int d,dim,i,maxident,nident,m,i1,i2,j,r,ll,ct,warn; molssptr mols; moleculeptr mptr; wallptr *wlist; char **name; boxssptr boxs; rxnptr rxn; double minboxsize,vol,amax,vol2,vol3; printf("PARAMETER CHECK\n"); warn=0; dim=sim->dim; maxident=sim->maxident; nident=sim->nident; mols=sim->mols; wlist=sim->wlist; name=sim->name; boxs=sim->boxs; for(ll=0;ll<2;ll++) for(m=0;mnl[ll];m++) { // molecules outside system mptr=mols->live[ll][m]; for(d=0;dpos[d]pos||mptr->pos[d]>wlist[2*d+1]->pos) { printf(" WARNING: molecule of type '%s' is outside system volume\n",name[mptr->ident]); warn++; }} minboxsize=boxs->size[0]; // too large rms steps for(d=1;dsize[d]size[d]; for(i=1;idifstep[i]>minboxsize/2.0) { printf(" WARNING: rms step length of %s is large compared to virtual box dimensions\n",name[i]); printf(" Try either more molecules per box or a shorter time step.\n"); warn++; } rxn=sim->rxn[2]; // difm ignored for reaction rates if(rxn) for(i1=1;i1difm[i1]) for(i2=1;i2nrxn[i=i1*maxident+i2];j++) if(rxn->rate[rxn->table[i][j]]) { printf(" WARNING: diffusion matrix for %s was ignored for calculating rate for reaction %s\n",name[i1],rxn->rname[rxn->table[i][j]]); warn++; } for(d=0;dtype=='p'&&wlist[2*d+1]->type!='p') { printf(" WARNING: only one wall on dimension %i has a periodic boundary condition\n",d); warn++; } rxn=sim->rxn[0]; // order 0 reactions if(rxn) for(r=0;rtotal;r++) if(rxn->rate2[r]<0&&rxn->rname[r][0]) { printf(" WARNING: reaction rate not set for reation order 0, name %s\n",rxn->rname[r]); warn++; } rxn=sim->rxn[1]; // order 1 reactions if(rxn) for(r=0;rtotal;r++) { if(rxn->rate2[r]<0&&rxn->rname[r][0]) { printf(" WARNING: reaction rate not set for reaction order 1, name %s\n",rxn->rname[r]); warn++; } else if(rxn->rate2[r]>0&&rxn->rate2[r]<10.0/RAND_MAX_30) { printf(" WARNING: order 1 reaction %s probability is at lower end of random number generator resolution\n",rxn->rname[r]); warn++; } else if(rxn->rate2[r]>(RAND_MAX_30-10.0)/RAND_MAX_30&&rxn->rate2[r]<1.0) { printf(" WARNING: order 1 reaction %s probability is at upper end of random number generator resolution\n",rxn->rname[r]); warn++; } else if(rxn->rate2[r]>0.2) { printf(" WARNING: order 1 reaction %s probability is quite high\n",rxn->rname[r]); warn++; }} rxn=sim->rxn[2]; // order 2 reactions if(rxn) { for(r=0;rtotal;r++) { if(rxn->rate2[r]<0&&rxn->rname[r][0]) { printf(" WARNING: reaction rate not set for reaction order 2, name %s\n",rxn->rname[r]); warn++; } else if(sqrt(rxn->rate2[r])>minboxsize) { printf(" WARNING: binding radius for order 2 reaction %s is larger than box dimensions\n",rxn->rname[r]); warn++; }} vol=1.0; for(d=0;dpos-wlist[2*d]->pos; vol2=0; for(i=1;inrxn[i*maxident+i1];j++) { r=rxn->table[i*maxident+i1][j]; if(amaxrate2[r])) amax=sqrt(rxn->rate2[r]); } ll=(mols->difc[i]==0); ct=0; for(m=0;mnl[ll];m++) if(mols->live[ll][m]->ident==i) ct++; vol3=ct*4.0/3.0*PI*amax*amax*amax; vol2+=vol3; if(vol3>vol/10.0) { printf(" WARNING: reactive volume of %s is %lf %% of total volume\n",sim->name[i],vol3/vol*100); warn++; }} if(vol2>vol/10.0) { printf(" WARNING: total reactive volume is a large fraction of total volume\n"); warn++; }} if(warn>=10) printf(" %i total warnings\n",warn); else if(!warn) printf(" No warnings\n"); printf("\n"); return; }