1 /* Steven Andrews, started 10/22/2001.
2 This is a library of functions for the Smoldyn program.
3 See documentation called SmoldynManual.pdf and SmoldynCodeDoc.pdf, and the Smoldyn
4 website, which is at www.smoldyn.org.
5 Copyright 2003-2016 by Steven Andrews. This work is distributed under the terms
6 of the Gnu Lesser General Public License (LGPL). */
7
8 #include <cfloat>
9 #include <cstdlib>
10 #include <cmath>
11 #include <cstdarg>
12 #include <cstdio>
13 #include <cstring>
14 #include <ctime>
15
16 #include "List.h"
17 #include "math2.h"
18 #include "opengl2.h"
19 #include "random2.h"
20 #include "Rn.h"
21 #include "SimCommand.h"
22 #include "string2.h"
23 #include "Zn.h"
24 #include <string>
25 #include <sstream>
26
27 #include "smoldyn.h"
28 #include "smoldynfuncs.h"
29 #include "smoldynconfigure.h"
30
31 #ifdef OPTION_VCELL
32 using std::stringstream;
33 #endif
34
35 #ifdef ENABLE_PYTHON_CALLBACK
36 #include "../python/CallbackFunc.h"
37 #endif
38
39 /******************************************************************************/
40 /************************* Global variable definitions ************************/
41 /******************************************************************************/
42
43 void (*LoggingCallback)(simptr,int,const char*,...)=NULL;
44 int ThrowThreshold=11;
45 FILE *LogFile=NULL;
46
47 //
48 // ErrorLineAndString contains ErrorString (a few places sprinf is used to put
49 // ErrorLine inside ErrorLineAndString therefore the size of ErrorLineAndString
50 // has to be bigger than ErrorString else compiler emits warning (for a good
51 // reason).
52 //
53 char ErrorLineAndString[STRCHARLONG]="";
54 char ErrorString[STRCHARLONG-STRCHAR-100]="";
55
56 int ErrorType=0;
57 char SimFlags[STRCHAR]="";
58 int VCellDefined=0;
59
60
61 /******************************************************************************/
62 /***************************** Simulation structure ***************************/
63 /******************************************************************************/
64
65
66 /******************************************************************************/
67 /****************************** Local declarations ****************************/
68 /******************************************************************************/
69
70 // error handling
71
72 // enumerated types
73 char *simss2string(enum SmolStruct ss,char *string);
74
75 // low level utilities
76
77 // memory management
78
79 // data structure output
80
81 // structure set up
82
83 // core simulation functions
84
85
86 /******************************************************************************/
87 /********************************* Error handling *****************************/
88 /******************************************************************************/
89
90
91 /* simSetLogging */
simSetLogging(FILE * logfile,void (* logFunction)(simptr,int,const char *,...))92 void simSetLogging(FILE *logfile,void (*logFunction)(simptr,int,const char*,...)) {
93 LogFile=logfile;
94 LoggingCallback=logFunction;
95 return; }
96
97
98 /* simSetThrowing */
simSetThrowing(int corethreshold)99 void simSetThrowing(int corethreshold) {
100 ThrowThreshold=corethreshold;
101 return; }
102
103
104 /* simLog */
simLog(simptr sim,int importance,const char * format,...)105 void simLog(simptr sim,int importance,const char* format, ...) {
106 char message[STRCHARLONG],*flags;
107 va_list arguments;
108 int qflag,vflag,wflag,sflag;
109 FILE *fptr;
110
111 va_start(arguments, format);
112 vsprintf(message, format, arguments);
113 va_end(arguments);
114
115 if(LoggingCallback)
116 (*LoggingCallback)(sim,importance,message);
117
118 if(sim && sim->logfile) fptr=sim->logfile;
119 else if(LogFile) fptr=LogFile;
120 else fptr=stdout;
121 if(sim) flags=sim->flags;
122 else flags=SimFlags;
123 qflag=strchr(flags,'q')?1:0;
124 sflag=strchr(flags,'s')?1:0;
125 vflag=strchr(flags,'v')?1:0;
126 wflag=strchr(flags,'w')?1:0;
127 if(!sflag)
128 if(vflag || importance>=2)
129 if(!qflag || importance>=5)
130 if(!wflag || importance<=4 || importance>=7)
131 fprintf(fptr,"%s", message);
132
133 //?? if(importance>=ThrowThreshold) throw message; //?? This is removed for now because I can't statically link Libsmoldyn if it's enabled
134
135 return; }
136
137
138 /* simParseError */
simParseError(simptr sim,ParseFilePtr pfp)139 void simParseError(simptr sim,ParseFilePtr pfp) {
140 char parseerrstr[STRCHAR],matherrstr[STRCHAR];
141
142 if(pfp) {
143 Parse_ReadFailure(pfp,parseerrstr);
144 sprintf(ErrorLineAndString,"%s\nMessage: %s",parseerrstr,ErrorString);
145 simLog(sim,8,"%s\nMessage: %s\n",parseerrstr,ErrorString);
146 if(strmatherror(matherrstr,1))
147 simLog(sim,8,"math error: %s\n",matherrstr); }
148 else {
149 sprintf(ErrorLineAndString,"%s",ErrorString);
150 simLog(sim,8,"%s",ErrorString); }
151 return; }
152
153
154 /******************************************************************************/
155 /********************************* enumerated types ***************************/
156 /******************************************************************************/
157
158
159 /* simstring2ss */
simstring2ss(char * string)160 enum SmolStruct simstring2ss(char *string) {
161 enum SmolStruct ans;
162
163 if(!strcmp(string,"molecule")) ans=SSmolec;
164 else if(!strcmp(string,"wall")) ans=SSwall;
165 else if(!strcmp(string,"reaction")) ans=SSrxn;
166 else if(!strcmp(string,"surface")) ans=SSsurf;
167 else if(!strcmp(string,"box")) ans=SSbox;
168 else if(!strcmp(string,"compartment")) ans=SScmpt;
169 else if(!strcmp(string,"port")) ans=SSport;
170 else if(!strcmp(string,"filament")) ans=SSfilament;
171 else if(!strcmp(string,"command")) ans=SScmd;
172 else if(!strcmp(string,"simulation")) ans=SSsim;
173 else if(!strcmp(string,"check")) ans=SScheck;
174 else if(!strcmp(string,"all")) ans=SSall;
175 else ans=SSnone;
176 return ans; }
177
178
179 /* simss2string */
simss2string(enum SmolStruct ss,char * string)180 char *simss2string(enum SmolStruct ss,char *string) {
181 if(ss==SSmolec) strcpy(string,"molecule");
182 else if(ss==SSwall) strcpy(string,"wall");
183 else if(ss==SSrxn) strcpy(string,"reaction");
184 else if(ss==SSsurf) strcpy(string,"surface");
185 else if(ss==SSbox) strcpy(string,"box");
186 else if(ss==SScmpt) strcpy(string,"compartment");
187 else if(ss==SSport) strcpy(string,"port");
188 else if(ss==SSfilament) strcpy(string,"filament");
189 else if(ss==SScmd) strcpy(string,"command");
190 else if(ss==SSsim) strcpy(string,"simulation");
191 else if(ss==SScheck) strcpy(string,"check");
192 else if(ss==SSall) strcpy(string,"all");
193 else strcpy(string,"none");
194 return string; }
195
196
197 /* simsc2string */
simsc2string(enum StructCond sc,char * string)198 char *simsc2string(enum StructCond sc,char *string) {
199 if(sc==SCinit) strcpy(string,"not initialized");
200 else if(sc==SClists) strcpy(string,"lists need updating");
201 else if(sc==SCparams) strcpy(string,"parameters need updating");
202 else if(sc==SCok) strcpy(string,"fully updated");
203 else strcpy(string,"none");
204 return string; }
205
206
207 /******************************************************************************/
208 /****************************** low level utilities ***************************/
209 /******************************************************************************/
210
211
212 /* simversionnumber */
simversionnumber(void)213 double simversionnumber(void) {
214 double v;
215
216 v=0;
217 sscanf(VERSION,"%lg",&v);
218 return v; }
219
220
221 /* Simsetrandseed */
Simsetrandseed(simptr sim,long int randseed)222 void Simsetrandseed(simptr sim,long int randseed) {
223 if(!sim) return;
224 sim->randseed=randomize(randseed);
225 return; }
226
227
228 /******************************************************************************/
229 /******************************* memory management ****************************/
230 /******************************************************************************/
231
232 /* simalloc */
simalloc(const char * fileroot)233 simptr simalloc(const char *fileroot) {
234 simptr sim;
235 int order;
236 enum EventType et;
237
238 sim=NULL;
239 CHECKMEM(sim=(simptr) malloc(sizeof(struct simstruct)));
240 sim->logfile=NULL;
241 sim->condition=SCinit;
242 sim->filepath=NULL;
243 sim->filename=NULL;
244 sim->flags=NULL;
245 sim->clockstt=time(NULL);
246 sim->elapsedtime=0;
247 Simsetrandseed(sim,-1);
248 for(et=(EventType)0;et<ETMAX;et=(EventType)(et+1)) sim->eventcount[et]=0;
249 sim->maxvar=0;
250 sim->nvar=0;
251 sim->varnames=NULL;
252 sim->varvalues=NULL;
253 sim->dim=0;
254 sim->accur=10;
255 sim->time=0;
256 sim->tmin=0;
257 sim->tmax=DBL_MAX;
258 sim->tbreak=DBL_MAX;
259 sim->dt=0;
260 sim->quitatend=0;
261 for(order=0;order<MAXORDER;order++) sim->rxnss[order]=NULL;
262 sim->ruless=NULL;
263 sim->mols=NULL;
264 sim->wlist=NULL;
265 sim->srfss=NULL;
266 sim->boxs=NULL;
267 sim->cmptss=NULL;
268 sim->portss=NULL;
269 sim->latticess=NULL;
270 sim->bngss=NULL;
271 sim->filss=NULL;
272 sim->cmds=NULL;
273 sim->graphss=NULL;
274
275 sim->diffusefn=&diffuse;
276 sim->surfaceboundfn=&checksurfacebound;
277 sim->surfacecollisionsfn=&checksurfaces;
278 sim->assignmols2boxesfn=&reassignmolecs;
279 sim->zeroreactfn=&zeroreact;
280 sim->unimolreactfn=&unireact;
281 sim->bimolreactfn=&bireact;
282 sim->checkwallsfn=&checkwalls;
283
284 CHECKMEM(sim->filepath=EmptyString());
285 CHECKMEM(sim->filename=EmptyString());
286 CHECKMEM(sim->flags=EmptyString());
287 CHECKMEM(sim->cmds=scmdssalloc(&docommand,(void*)sim,fileroot));
288
289 simsetvariable(sim,"time",sim->time);
290 simsetvariable(sim,"x",dblnan());
291 simsetvariable(sim,"y",dblnan());
292 simsetvariable(sim,"z",dblnan());
293 simsetvariable(sim,"r",dblnan());
294
295 #if ENABLE_PYTHON_CALLBACK
296 sim->ncallbacks = 0;
297 sim->simstep = 0;
298 #endif
299
300 return sim;
301
302 failure:
303 simfree(sim);
304 simLog(NULL,10,"Unable to allocate memory in simalloc");
305 return NULL; }
306
307
308 /* simfree */
simfree(simptr sim)309 void simfree(simptr sim) {
310 int dim,order,maxsrf;
311
312 if(!sim) return;
313 dim=sim->dim;
314 maxsrf=sim->srfss?sim->srfss->maxsrf:0;
315
316 graphssfree(sim->graphss);
317 scmdssfree(sim->cmds);
318 filssfree(sim->filss);
319 latticessfree(sim->latticess);
320 portssfree(sim->portss);
321 compartssfree(sim->cmptss);
322 boxssfree(sim->boxs);
323 surfacessfree(sim->srfss);
324 wallsfree(sim->wlist,dim);
325 molssfree(sim->mols,maxsrf);
326 rulessfree(sim->ruless);
327 for(order=0;order<MAXORDER;order++) rxnssfree(sim->rxnss[order]);
328
329 for(size_t v=0;v<(size_t)sim->maxvar;v++)
330 free(sim->varnames[v]);
331 free(sim->varnames);
332
333 #ifdef ENABLE_PYTHON_CALLBACK
334 for(size_t v = 0; v < sim->ncallbacks; v++)
335 if(sim->callbacks[v])
336 free(sim->callbacks[v]);
337 #endif
338
339 free(sim->varvalues);
340 free(sim->flags);
341 free(sim->filename);
342 free(sim->filepath);
343 free(sim);
344 return; }
345
346
347 /* simfuncfree */
simfuncfree(void)348 void simfuncfree(void) {
349 strEnhWildcardMatch(NULL,NULL);
350 strEnhWildcardMatchAndSub(NULL,NULL,NULL,NULL);
351 strevalfunction(NULL,NULL,NULL,NULL,NULL,NULL,0);
352 return; }
353
354
355 /* sim_expandvariables */
simexpandvariables(simptr sim,int spaces)356 int simexpandvariables(simptr sim,int spaces) {
357 char **newvarnames;
358 double *newvarvalues;
359 int newmaxvar,i,newnvar;
360
361 newmaxvar=sim->maxvar+spaces;
362 newvarnames=(char**) calloc(newmaxvar,sizeof (char*));
363 if(!newvarnames) return 1;
364 newvarvalues=(double*) calloc(newmaxvar,sizeof (double));
365 if(!newvarvalues) return 1;
366 for(i=0;i<sim->nvar && i<newmaxvar;i++) newvarnames[i]=sim->varnames[i];
367 for(i=0;i<sim->nvar && i<newmaxvar;i++) newvarvalues[i]=sim->varvalues[i];
368 newnvar=i;
369 for(i=newnvar;i<newmaxvar;i++) {
370 newvarnames[i]=EmptyString();
371 if(!newvarnames[i]) return 1; }
372 for(i=newnvar;i<newmaxvar;i++) newvarvalues[i]=0;
373 free(sim->varnames);
374 free(sim->varvalues);
375 sim->varnames=newvarnames;
376 sim->varvalues=newvarvalues;
377 sim->maxvar=newmaxvar;
378 sim->nvar=newnvar;
379 return 0; }
380
381
382 /******************************************************************************/
383 /***************************** data structure output **************************/
384 /******************************************************************************/
385
386
387 /* simoutput */
simoutput(simptr sim)388 void simoutput(simptr sim) {
389 int v;
390
391 simLog(sim,2,"SIMULATION PARAMETERS\n");
392 if(!sim) {
393 simLog(sim,2," No simulation parameters\n\n");
394 return; }
395 if(sim->filename[0]!='\0')
396 simLog(sim,2," file: %s%s\n",sim->filepath,sim->filename);
397 simLog(sim,2," starting clock time: %s",ctime(&sim->clockstt));
398 simLog(sim,2," %i dimensions\n",sim->dim);
399 if(sim->accur<10) simLog(sim,2," Accuracy level: %g\n",sim->accur);
400 else simLog(sim,1," Accuracy level: %g\n",sim->accur);
401 simLog(sim,2," Random number seed: %li\n",sim->randseed);
402 simLog(sim,sim->nvar>5?2:1," %i variable%s defined:\n",sim->nvar,sim->nvar>1?"s":"");
403 for(v=0;v<sim->nvar && v<5;v++)
404 simLog(sim,1," %s = %g\n",sim->varnames[v],sim->varvalues[v]);
405 for(;v<sim->nvar;v++)
406 simLog(sim,2," %s = %g\n",sim->varnames[v],sim->varvalues[v]);
407
408 simLog(sim,2," Time from %g to %g step %g\n",sim->tmin,sim->tmax,sim->dt);
409 if(sim->time!=sim->tmin) simLog(sim,2," Current time: %g\n",sim->time);
410 simLog(sim,2,"\n");
411 return; }
412
413
414 /* simsystemoutput */
simsystemoutput(simptr sim)415 void simsystemoutput(simptr sim) {
416 int order,vflag;
417
418 if(!sim) {
419 simLog(sim,2," No simulation parameters\n\n");
420 return; }
421 vflag=strchr(sim->flags,'v')?1:0;
422 simoutput(sim);
423 graphssoutput(sim);
424 walloutput(sim);
425 molssoutput(sim);
426 surfaceoutput(sim);
427 scmdoutput(sim->cmds);
428 boxssoutput(sim);
429 if(vflag)
430 // boxoutput(sim->boxs,0,sim->boxs->nbox,sim->dim);
431 boxoutput(sim->boxs,0,10,sim->dim);
432 for(order=0;order<MAXORDER;order++)
433 rxnoutput(sim,order);
434 ruleoutput(sim);
435 compartoutput(sim);
436 portoutput(sim);
437 bngoutput(sim);
438 latticeoutput(sim);
439 filssoutput(sim);
440 return; }
441
442
443 /* writesim */
writesim(simptr sim,FILE * fptr)444 void writesim(simptr sim,FILE *fptr) {
445 fprintf(fptr,"# General simulation parameters\n");
446 fprintf(fptr,"# Configuration file: %s%s\n",sim->filepath,sim->filename);
447 fprintf(fptr,"dim %i\n",sim->dim);
448 fprintf(fptr,"# random_seed for prior simulation was %li\n",sim->randseed);
449 fprintf(fptr,"random_seed %li # this is a new random number\n",(long int)randULI());
450 fprintf(fptr,"time_start %g\n",sim->tmin);
451 fprintf(fptr,"time_stop %g\n",sim->tmax);
452 fprintf(fptr,"time_step %g\n",sim->dt);
453 fprintf(fptr,"time_now %g\n",sim->time);
454 fprintf(fptr,"accuracy %g\n",sim->accur);
455 if(sim->boxs->mpbox) fprintf(fptr,"molperbox %g\n",sim->boxs->mpbox);
456 else if(sim->boxs->boxsize) fprintf(fptr,"boxsize %g\n",sim->boxs->boxsize);
457 fprintf(fptr,"\n");
458 return; }
459
460
461 /* checksimparams */
checksimparams(simptr sim)462 int checksimparams(simptr sim) {
463 int warn,error,warndiff;
464 char string[STRCHAR];
465
466 simLog(sim,2,"PARAMETER CHECK\n");
467 warn=error=0;
468
469 error+=checkmolparams(sim,&warndiff);warn+=warndiff;
470 error+=checkboxparams(sim,&warndiff);warn+=warndiff;
471 error+=checkwallparams(sim,&warndiff);warn+=warndiff;
472 error+=checkrxnparams(sim,&warndiff);warn+=warndiff;
473 error+=checkruleparams(sim,&warndiff);warn+=warndiff;
474 error+=checksurfaceparams(sim,&warndiff);warn+=warndiff;
475 error+=checkcompartparams(sim,&warndiff);warn+=warndiff;
476 error+=checkportparams(sim,&warndiff);warn+=warndiff;
477 error+=checklatticeparams(sim,&warndiff);warn+=warndiff;
478 error+=filcheckparams(sim,&warndiff);warn+=warndiff;
479 error+=checkgraphicsparams(sim,&warndiff);warn+=warndiff;
480 error+=checkbngparams(sim,&warndiff);warn+=warndiff;
481
482 if(sim->condition!=SCok) {
483 warn++;
484 simLog(sim,7," WARNING: simulation structure %s\n",simsc2string(sim->condition,string)); }
485
486 if(error>0) simLog(sim,2," %i total errors\n",error);
487 else simLog(sim,2," No errors\n");
488 if(warn>0) simLog(sim,2," %i total warnings\n",warn);
489 else simLog(sim,2," No warnings\n");
490 simLog(sim,2,"\n");
491 return error; }
492
493
494
495 /******************************************************************************/
496 /******************************* structure set up *****************************/
497 /******************************************************************************/
498
499
500 /* simsetcondition */
simsetcondition(simptr sim,enum StructCond cond,int upgrade)501 void simsetcondition(simptr sim,enum StructCond cond,int upgrade) {
502 if(!sim) return;
503 if(upgrade==0 && sim->condition>cond) sim->condition=cond;
504 else if(upgrade==1 && sim->condition<cond) sim->condition=cond;
505 else if(upgrade==2) sim->condition=cond;
506 return; }
507
508
509 /* simsetvariable */
simsetvariable(simptr sim,const char * name,double value)510 int simsetvariable(simptr sim,const char *name,double value) {
511 int v,er;
512
513 v=stringfind(sim->varnames,sim->nvar,name);
514 if(v<0) {
515 if(sim->nvar==sim->maxvar) {
516 er=simexpandvariables(sim,2+2*sim->nvar);
517 if(er) return er; }
518 v=sim->nvar++;
519 strcpy(sim->varnames[v],name); }
520 sim->varvalues[v]=value;
521 return 0; }
522
523
524 /* simsetdim */
simsetdim(simptr sim,int dim)525 int simsetdim(simptr sim,int dim) {
526 if(sim->dim!=0) return 2;
527 if(dim<1 || dim>DIMMAX) return 3;
528 sim->dim=dim;
529 return 0; }
530
531
532 /* simsettime */
simsettime(simptr sim,double time,int code)533 int simsettime(simptr sim,double time,int code) {
534 int er;
535 static int timedefined=0;
536
537 if(code==-1) return timedefined;
538 else if(code==0) timedefined|=1;
539 else if(code==1) timedefined|=2;
540 else if(code==2) timedefined|=4;
541 else if(code==3) timedefined|=8;
542 else if(code==4) timedefined|=16;
543
544 er=0;
545 if(code==0) {
546 sim->time=time;
547 simsetvariable(sim,"time",time); }
548 else if(code==1) sim->tmin=time;
549 else if(code==2) sim->tmax=time;
550 else if(code==3) {
551 if(time>0) {
552 sim->dt=time;
553 molsetcondition(sim->mols,SCparams,0);
554 rxnsetcondition(sim,-1,SCparams,0);
555 surfsetcondition(sim->srfss,SCparams,0);
556 scmdsetcondition(sim->cmds,1,0); }
557 else er=2; }
558 else if(code==4) sim->tbreak=time;
559 else er=1;
560 return er; }
561
562
563 /* simreadstring */
simreadstring(simptr sim,ParseFilePtr pfp,const char * word,char * line2)564 int simreadstring(simptr sim,ParseFilePtr pfp,const char *word,char *line2) {
565 char nm[STRCHAR],nm1[STRCHAR],shapenm[STRCHAR],ch,rname[STRCHAR],fname[STRCHAR],pattern[STRCHAR];
566 char str1[STRCHAR],str2[STRCHAR],str3[STRCHAR];
567 char errstr[STRCHARLONG];
568 int er,i,nmol,d,i1,s,c,ll,order,*index,ft;
569 int rulelist[MAXORDER+MAXPRODUCT],r,ord,rct,prd,itct,prt,lt,detailsi[8];
570 long int pserno,sernolist[MAXPRODUCT];
571 double flt1,flt2,v1[DIMMAX*DIMMAX],v2[4],poslo[DIMMAX],poshi[DIMMAX],thick;
572 enum MolecState ms,rctstate[MAXORDER];
573 enum PanelShape ps;
574 enum RevParam rpart;
575 enum LightParam ltparam;
576 enum SpeciesRepresentation replist[MAXORDER+MAXPRODUCT];
577 rxnptr rxn;
578 compartptr cmpt;
579 surfaceptr srf;
580 portptr port;
581 filamentptr fil;
582 filamenttypeptr filtype;
583 filamentssptr filss;
584 long int li1;
585 listptrli lilist;
586 listptrv vlist;
587
588 int dim,nvar;
589 char **varnames;
590 double *varvalues;
591
592 dim=sim->dim;
593 nvar=sim->nvar;
594 varnames=sim->varnames;
595 varvalues=sim->varvalues;
596
597 // variables
598
599 if(!strcmp(word,"variable")) { // variable
600 itct=sscanf(line2,"%s",nm);
601 CHECKS(itct==1,"variable format: name = value");
602 CHECKS(strcmp(nm,"time"),"'time' cannot be used as a variable name; it is pre-defined as the simulation time");
603 CHECKS(strcmp(nm,"x") && strcmp(nm,"y") && strcmp(nm,"z") && strcmp(nm,"r"),"x, y, z, and r are reserved variable names");
604 CHECKS(strokname(nm),"variable name has to start with a letter and then be alphanumeric with optional underscores");
605 line2=strnword(line2,2);
606 CHECKS(line2,"variable format: name = value");
607 if(line2[0]=='=') {
608 line2=strnword(line2,2);
609 CHECKS(line2,"variable format: name = value"); }
610 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
611 CHECKS(itct==1,"cannot read variable value");
612 er=simsetvariable(sim,nm,flt1);
613 CHECKS(!er,"out of memory allocating variable space");
614 CHECKS(!strnword(line2,2),"unexpected text following variable"); }
615
616 // space and time
617
618 else if(!strcmp(word,"dim")) { // dim
619 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&d);
620 CHECKS(itct==1,"dim needs to be an integer");
621 er=simsetdim(sim,d);
622 CHECKS(er!=2,"dim can only be entered once");
623 CHECKS(er!=3,"dim has to be between 1 and 3");
624 dim=sim->dim;
625 CHECKS(!strnword(line2,2),"unexpected text following dim"); }
626
627 else if(!strcmp(word,"boundaries")) { // boundaries
628 CHECKS(dim>0,"need to enter dim before boundaries");
629 itct=strmathsscanf(line2,"%s %mlg %mlg",varnames,varvalues,nvar,nm,&flt1,&flt2);
630 CHECKS(itct==3,"boundaries format: dimension position position [type]");
631 if(!strcmp(nm,"0") || !strcmp(nm,"x")) d=0;
632 else if(!strcmp(nm,"1") || !strcmp(nm,"y")) d=1;
633 else if(!strcmp(nm,"2") || !strcmp(nm,"z")) d=2;
634 else d=3;
635 CHECKS(d>=0 && d<dim,"illegal boundaries dimension value");
636 CHECKS(flt1<flt2,"the first boundary position needs to be smaller than the second one");
637 line2=strnword(line2,4);
638 ch='t';
639 if(line2) {
640 itct=sscanf(line2,"%c",&ch);
641 CHECKS(itct==1,"boundary type character could not be read");
642 CHECKS(ch=='r' || ch=='p' || ch=='a' || ch=='t',"boundaries type needs to be r, p, a, or t");
643 line2=strnword(line2,2); }
644 er=walladd(sim,d,0,flt1,ch);
645 CHECKS(er!=1,"out of memory adding wall");
646 CHECKS(er!=2,"SMOLDYN BUG: adding wall");
647 er=walladd(sim,d,1,flt2,ch);
648 CHECKS(er!=1,"out of memory adding wall");
649 CHECKS(er!=2,"SMOLDYN BUG: adding wall");
650 CHECKS(!line2,"unexpected text following boundaries"); }
651
652 else if(!strcmp(word,"low_wall")) { // low_wall
653 CHECKS(dim>0,"need to enter dim before low_wall");
654 itct=strmathsscanf(line2,"%s %mlg %c",varnames,varvalues,nvar,nm,&flt1,&ch);
655 CHECKS(itct==3,"low_wall format: dimension position type");
656 if(!strcmp(nm,"0") || !strcmp(nm,"x")) d=0;
657 else if(!strcmp(nm,"1") || !strcmp(nm,"y")) d=1;
658 else if(!strcmp(nm,"2") || !strcmp(nm,"z")) d=2;
659 else d=3;
660 CHECKS(d>=0 && d<dim,"low_wall dimension needs to be between 0 and dim-1");
661 CHECKS(ch=='r' || ch=='p' || ch=='a' || ch=='t',"low_wall type needs to be r, p, a, or t");
662 er=walladd(sim,d,0,flt1,ch);
663 CHECKS(er!=1,"out of memory adding wall");
664 CHECKS(er!=2,"SMOLDYN BUG: adding wall");
665 CHECKS(!strnword(line2,4),"unexpected text following low_wall"); }
666
667 else if(!strcmp(word,"high_wall")) { // high_wall
668 CHECKS(dim>0,"need to enter dim before high_wall");
669 itct=strmathsscanf(line2,"%s %mlg %c",varnames,varvalues,nvar,nm,&flt1,&ch);
670 CHECKS(itct==3,"high_wall format: dimension position type");
671 if(!strcmp(nm,"0") || !strcmp(nm,"x")) d=0;
672 else if(!strcmp(nm,"1") || !strcmp(nm,"y")) d=1;
673 else if(!strcmp(nm,"2") || !strcmp(nm,"z")) d=2;
674 else d=3;
675 CHECKS(d>=0 && d<dim,"high_wall dimension needs to be between 0 and dim-1");
676 CHECKS(ch=='r' || ch=='p' || ch=='a' || ch=='t',"high_wall type needs to be r, p, a, or t");
677 er=walladd(sim,d,1,flt1,ch);
678 CHECKS(er!=1,"out of memory adding wall");
679 CHECKS(er!=2,"SMOLDYN BUG: adding wall");
680 CHECKS(!strnword(line2,4),"unexpected text following high_wall"); }
681
682 else if(!strcmp(word,"time_start")) { // time_start
683 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
684 CHECKS(itct==1,"time_start needs to be a number");
685 simsettime(sim,flt1,0);
686 simsettime(sim,flt1,1);
687 CHECKS(!strnword(line2,2),"unexpected text following time_start"); }
688
689 else if(!strcmp(word,"time_stop")) { // time_stop
690 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
691 CHECKS(itct==1,"time_stop needs to be a number");
692 simsettime(sim,flt1,2);
693 CHECKS(!strnword(line2,2),"unexpected text following time_stop"); }
694
695 else if(!strcmp(word,"time_step")) { // time_step
696 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
697 CHECKS(itct==1,"time_step needs to be a number");
698 er=simsettime(sim,flt1,3);
699 CHECKS(!er,"time step must be >0");
700 CHECKS(!strnword(line2,2),"unexpected text following time_step"); }
701
702 else if(!strcmp(word,"time_now")) { // time_now
703 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
704 CHECKS(itct==1,"time_now needs to be a number");
705 simsettime(sim,flt1,0);
706 CHECKS(!strnword(line2,2),"unexpected text following time_now"); }
707
708 // molecules
709
710 else if(!strcmp(word,"max_species") || !strcmp(word,"max_names")) { // max_species, max_names
711 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&i1);
712 CHECKS(itct==1,"max_species needs to be an integer");
713 er=molenablemols(sim,i1);
714 CHECKS(er!=1,"out of memory");
715 CHECKS(er!=2,"cannot decrease the number of allocated species");
716 CHECKS(!strnword(line2,2),"unexpected text following max_species"); }
717
718 else if(!strcmp(word,"species") || !strcmp(word,"names") || !strcmp(word,"name")) {// species, names, name
719 while(line2) {
720 itct=sscanf(line2,"%s",nm);
721 CHECKS(itct==1,"failed to read species name");
722 er=moladdspecies(sim,nm);
723 CHECKS(er!=-1,"failed to allocate memory");
724 CHECKS(er!=-4,"'empty' is not a permitted species name");
725 CHECKS(er!=-5,"this species has already been declared");
726 CHECKS(er!=-6,"'?' and '*' are not permitted in species names");
727 line2=strnword(line2,2); }}
728
729 else if(!strcmp(word,"species_group")) { // species_group
730 itct=sscanf(line2,"%s",nm);
731 CHECKS(itct==1,"species_group format: group species species ...");
732 line2=strnword(line2,2);
733 if(!line2) {
734 er=moladdspeciesgroup(sim,nm,NULL,0);
735 CHECKS(er!=-4,"species name not recognized");
736 CHECKS(er!=-6,"cannot read wildcard logic");
737 CHECKS(er!=-7,"error allocating memory");
738 CHECKS(er!=-8,"molecule states are not permitted");
739 CHECKS(er!=-9,"species group name cannot be a species name");
740 CHECKS(!er,"BUG: failed to add species group"); }
741 else
742 while(line2) {
743 itct=sscanf(line2,"%s",nm1);
744 CHECKS(itct==1,"species_group format: group species species ...");
745 er=moladdspeciesgroup(sim,nm,nm1,0);
746 CHECKS(er!=-4,"species name not recognized");
747 CHECKS(er!=-6,"cannot read wildcard logic");
748 CHECKS(er!=-7,"error allocating memory");
749 CHECKS(er!=-8,"molecule states are not permitted");
750 CHECKS(!er,"BUG: failed to add species group");
751 line2=strnword(line2,2); }}
752
753 else if(!strcmp(word,"max_mol")) { // max_mol
754 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&i1);
755 CHECKS(itct==1,"max_mol needs to be an integer");
756 er=molsetmaxmol(sim,i1);
757 CHECKS(er!=1,"out of memory");
758 CHECKS(er!=5,"more molecule already exist than are requested with max_mol");
759 CHECKS(!strnword(line2,2),"unexpected text following max_mol"); }
760
761 else if(!strcmp(word,"difc")) { // difc
762 CHECKS(sim->mols,"need to enter species before difc");
763 er=molstring2index1(sim,line2,&ms,&index);
764 CHECKS(er!=-1,"difc format: species[(state)] value");
765 CHECKS(er!=-2,"mismatched or improper parentheses around molecule state");
766 CHECKS(er!=-3,"cannot read molecule state value");
767 CHECKS(er!=-4,"molecule name not recognized");
768 CHECKS(er!=-6,"cannot read wildcard logic");
769 CHECKS(er!=-7,"error allocating memory");
770 CHECKS(ms<MSMAX || ms==MSall,"invalid state");
771 CHECKS(line2=strnword(line2,2),"missing data for difc");
772 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
773 CHECKS(itct==1,"diffusion coefficient value cannot be read");
774 CHECKS(flt1>=0,"diffusion coefficient needs to be >=0");
775 molsetdifc(sim,0,index,ms,flt1);
776 CHECKS(!strnword(line2,2),"unexpected text following difc"); }
777
778 else if(!strcmp(word,"difc_rule")) { // difc_rule
779 er=molstring2pattern(line2,&ms,pattern,0);
780 CHECKS(er!=-1,"difc_rule format: species[(state)] value");
781 CHECKS(er!=-2,"mismatched or improper parentheses around molecule state");
782 CHECKS(er!=-3,"cannot read molecule state value");
783 CHECKS(ms<MSMAX || ms==MSall,"invalid state");
784 CHECKS(line2=strnword(line2,2),"missing data for difc_rule");
785 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
786 CHECKS(itct==1,"diffusion coefficient value cannot be read");
787 CHECKS(flt1>=0,"diffusion coefficient needs to be >=0");
788 er=RuleAddRule(sim,RTdifc,NULL,pattern,&ms,NULL,flt1,NULL,NULL);
789 CHECKS(!er,"out of memory in difc_rule");
790 CHECKS(!strnword(line2,2),"unexpected text following difc_rule"); }
791
792 else if(!strcmp(word,"difm")) { // difm
793 CHECKS(sim->mols,"need to enter species before difm");
794 er=molstring2index1(sim,line2,&ms,&index);
795 CHECKS(er!=-1,"difm format: name[(state)] values");
796 CHECKS(er!=-2,"mismatched or improper parentheses around molecule state");
797 CHECKS(er!=-3,"cannot read molecule state value");
798 CHECKS(er!=-4,"molecule name not recognized");
799 CHECKS(er!=-6,"cannot read wildcard logic");
800 CHECKS(er!=-7,"error allocating memory");
801 CHECKS(ms<MSMAX || ms==MSall,"invalid state");
802 CHECKS(line2=strnword(line2,2),"missing matrix in difm");
803 for(d=0;d<dim*dim;d++) {
804 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&v1[d]);
805 CHECKS(itct==1,"incomplete matrix in difm");
806 line2=strnword(line2,2); }
807 CHECKS(molsetdifm(sim,0,index,ms,v1)==0,"out of memory in difm");
808 CHECKS(!line2,"unexpected text following difm"); }
809
810 else if(!strcmp(word,"difm_rule")) { // difm_rule
811 er=molstring2pattern(line2,&ms,pattern,0);
812 CHECKS(er!=-1,"difm_rule format: name[(state)] values");
813 CHECKS(er!=-2,"mismatched or improper parentheses around molecule state");
814 CHECKS(er!=-3,"cannot read molecule state value");
815 CHECKS(ms<MSMAX || ms==MSall,"invalid state");
816 CHECKS(line2=strnword(line2,2),"missing matrix in difm");
817 for(d=0;d<dim*dim;d++) {
818 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&v1[d]);
819 CHECKS(itct==1,"incomplete matrix in difm");
820 line2=strnword(line2,2); }
821 er=RuleAddRule(sim,RTdifm,NULL,pattern,&ms,NULL,0,NULL,v1);
822 CHECKS(!er,"out of memory in difm_rule");
823 CHECKS(!line2,"unexpected text following difm_rule"); }
824
825 else if(!strcmp(word,"drift")) { // drift
826 CHECKS(sim->mols,"need to enter species before drift");
827 er=molstring2index1(sim,line2,&ms,&index);
828 CHECKS(er!=-1,"drift format: species[(state)] vector");
829 CHECKS(er!=-2,"mismatched or improper parentheses around molecule state");
830 CHECKS(er!=-3,"cannot read molecule state value");
831 CHECKS(er!=-4,"molecule name not recognized");
832 CHECKS(er!=-6,"cannot read wildcard logic");
833 CHECKS(er!=-7,"error allocating memory");
834 CHECKS(ms<MSMAX || ms==MSall,"invalid state");
835 CHECKS(line2=strnword(line2,2),"missing vector in drift");
836 for(d=0;d<dim;d++) {
837 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&v1[d]);
838 CHECKS(itct==1,"incomplete vector in drift");
839 line2=strnword(line2,2); }
840 CHECKS(molsetdrift(sim,0,index,ms,v1)==0,"out of memory in drift");
841 CHECKS(!line2,"unexpected text following drift"); }
842
843 else if(!strcmp(word,"drift_rule")) { // drift_rule
844 er=molstring2pattern(line2,&ms,pattern,0);
845 CHECKS(er!=-1,"drift_rule format: species[(state)] vector");
846 CHECKS(er!=-2,"mismatched or improper parentheses around molecule state");
847 CHECKS(er!=-3,"cannot read molecule state value");
848 CHECKS(ms<MSMAX || ms==MSall,"invalid state");
849 CHECKS(line2=strnword(line2,2),"missing vector in drift");
850 for(d=0;d<dim;d++) {
851 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&v1[d]);
852 CHECKS(itct==1,"incomplete vector in drift");
853 line2=strnword(line2,2); }
854 er=RuleAddRule(sim,RTdrift,NULL,pattern,&ms,NULL,0,NULL,v1);
855 CHECKS(!er,"out of memory in drift_rule");
856 CHECKS(!line2,"unexpected text following drift_rule"); }
857
858 else if(!strcmp(word,"surface_drift")) { // surface_drift
859 CHECKS(sim->dim>1,"system dimensionality needs to be 2 or 3");
860 CHECKS(sim->mols,"need to enter species before surface_drift");
861 CHECKS(sim->srfss,"need to enter surfaces before surface_drift");
862 CHECKS(sim->srfss->nsrf>0,"at least one surface needs to be defined before surface_drift");
863 er=molstring2index1(sim,line2,&ms,&index);
864 CHECKS(er!=-1,"surface_drift format: species(state) surface panel_shape values");
865 CHECKS(er!=-2,"mismatched or improper parentheses around molecule state");
866 CHECKS(er!=-3,"cannot read molecule state value");
867 CHECKS(er!=-4,"molecule name not recognized");
868 CHECKS(er!=-6,"cannot read wildcard logic");
869 CHECKS(er!=-7,"error allocating memory");
870 CHECKS(ms<MSMAX || ms==MSall,"invalid state");
871 CHECKS(ms!=MSsoln,"state needs to be surface-bound");
872 CHECKS(line2=strnword(line2,2),"surface_drift format: species(state) surface panel_shape values");
873 itct=sscanf(line2,"%s %s",nm,shapenm);
874 CHECKS(itct==2,"surface_drift format: species(state) surface panel_shape values");
875 if(!strcmp(nm,"all")) s=-1;
876 else {
877 s=stringfind(sim->srfss->snames,sim->srfss->nsrf,nm);
878 CHECKS(s>=0,"surface name in surface_drift is not recognized"); }
879 ps=surfstring2ps(shapenm);
880 CHECKS(ps!=PSnone,"in surface_drift, panel shape name not recognized");
881 CHECKS(line2=strnword(line2,3),"missing vector in surface_drift");
882 for(d=0;d<dim-1;d++) {
883 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&v1[d]);
884 CHECKS(itct==1,"incomplete vector in surface_drift");
885 line2=strnword(line2,2); }
886 CHECKS(molsetsurfdrift(sim,0,index,ms,s,ps,v1)==0,"out of memory in surface_drift");
887 CHECKS(!line2,"unexpected text following surface_drift"); }
888
889 else if(!strcmp(word,"surface_drift_rule")) { // surface_drift_rule
890 CHECKS(sim->dim>1,"system dimensionality needs to be 2 or 3");
891 CHECKS(sim->srfss,"need to enter surfaces before surface_drift_rule");
892 CHECKS(sim->srfss->nsrf>0,"at least one surface needs to be defined before surface_drift_rule");
893 er=molstring2pattern(line2,&ms,pattern,0);
894 CHECKS(er!=-1,"surface_drift_rule format: species(state) surface panel_shape values");
895 CHECKS(er!=-2,"mismatched or improper parentheses around molecule state");
896 CHECKS(er!=-3,"cannot read molecule state value");
897 CHECKS(ms<MSMAX || ms==MSall,"invalid state");
898 CHECKS(ms!=MSsoln,"state needs to be surface-bound");
899 CHECKS(line2=strnword(line2,2),"surface_drift_rule format: species(state) surface panel_shape values");
900 itct=sscanf(line2,"%s %s",nm,shapenm);
901 CHECKS(itct==2,"surface_drift_rule format: species(state) surface panel_shape values");
902 if(!strcmp(nm,"all")) s=-1;
903 else {
904 s=stringfind(sim->srfss->snames,sim->srfss->nsrf,nm);
905 CHECKS(s>=0,"surface name in surface_drift_rule is not recognized"); }
906 ps=surfstring2ps(shapenm);
907 CHECKS(ps!=PSnone,"in surface_drift_rule, panel shape name not recognized");
908 CHECKS(line2=strnword(line2,3),"missing vector in surface_drift");
909 for(d=0;d<dim-1;d++) {
910 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&v1[d]);
911 CHECKS(itct==1,"incomplete vector in surface_drift_rule");
912 line2=strnword(line2,2); }
913 detailsi[0]=s;
914 detailsi[1]=(int)ps;
915 er=RuleAddRule(sim,RTsurfdrift,NULL,pattern,&ms,NULL,0,detailsi,v1);
916 CHECKS(!er,"out of memory in surface_drift_rule");
917 CHECKS(!line2,"unexpected text following surface_drift_rule"); }
918
919 else if(!strcmp(word,"mol")) { // mol
920 CHECKS(sim->mols,"need to enter species before mol");
921 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&nmol);
922 CHECKS(itct==1,"mol format: number name position_vector");
923 CHECKS(nmol>=0,"number of molecules added needs to be >=0");
924 CHECKS(line2=strnword(line2,2),"mol format: number name position_vector");
925 i=molstring2index1(sim,line2,&ms,NULL);
926 CHECKS(i!=0,"wildcards and species groups are not permitted");
927 CHECKS(i!=-1,"error reading molecule name");
928 CHECKS(i!=-2,"mismatched or improper parentheses around molecule state");
929 CHECKS(i!=-3,"cannot read molecule state value");
930 CHECKS(i!=-4,"molecule name not recognized");
931 CHECKS(i!=-5,"'all' species is not permitted");
932 CHECKS(i!=-6,"cannot read wildcard logic");
933 CHECKS(i!=-7,"out of memory");
934 CHECKS(ms==MSsoln,"state needs to be solution");
935 if(sim->wlist)
936 systemcorners(sim,poslo,poshi);
937 else
938 for(d=0;d<dim;d++) {poslo[d]=0;poshi[d]=1;}
939 for(d=0;d<dim;d++) {
940 CHECKS(line2=strnword(line2,2),"insufficient position data for mol");
941 if(line2[0]=='u') {
942 CHECKS(sim->wlist,"need to enter boundaries before using 'u' in a mol statement"); }
943 else {
944 itct=sscanf(line2,"%lg-%lg",&flt1,&flt2);
945 if(itct==2) {
946 poslo[d]=flt1;
947 poshi[d]=flt2; }
948 else {
949 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
950 CHECKS(itct==1,"cannot read position value for mol");
951 poslo[d]=poshi[d]=flt1; }}}
952 er=addmol(sim,nmol,i,poslo,poshi,0);
953 CHECKS(!er,"more molecules assigned than permitted with max_mol");
954 CHECKS(!strnword(line2,2),"unexpected text following mol"); }
955
956 else if(!strcmp(word,"surface_mol")) { // surface_mol
957 CHECKS(sim->mols,"need to enter species before surface_mol");
958 CHECKS(sim->srfss,"surfaces need to be defined before surface_mol statement");
959 CHECKS(sim->srfss->nsrf>0,"at least one surface needs to be defined before surface_mol");
960 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&nmol);
961 CHECKS(itct==1,"surface_mol format: nmol species(state) surface panel_shape panel_name [position]");
962 CHECKS(nmol>=0,"in surface_mol, the number of molecules needs to be at least 0");
963 CHECKS(line2=strnword(line2,2),"surface_mol format: nmol species(state) surface panel_shape panel_name [position]");
964 i=molstring2index1(sim,line2,&ms,NULL);
965 CHECKS(i!=0,"wildcards and species groups are not permitted");
966 CHECKS(i!=-1,"error reading molecule name");
967 CHECKS(i!=-2,"mismatched or improper parentheses around molecule state");
968 CHECKS(i!=-3,"cannot read molecule state value");
969 CHECKS(i!=-4,"molecule name not recognized");
970 CHECKS(i!=-5,"'all' species is not permitted");
971 CHECKS(i!=-6,"cannot read wildcard logic");
972 CHECKS(i!=-7,"out of memory");
973 CHECKS(ms<MSMAX && ms!=MSsoln,"state needs to be front, back, up, or down");
974 CHECKS(line2=strnword(line2,2),"surface_mol format: nmol species(state) surface panel_shape panel_name [position]");
975 itct=sscanf(line2,"%s %s %s",nm,shapenm,nm1);
976 CHECKS(itct==3,"surface_mol format: nmol species(state) surface panel_shape panel_name [position]");
977 if(!strcmp(nm,"all")) s=-1;
978 else {
979 s=stringfind(sim->srfss->snames,sim->srfss->nsrf,nm);
980 CHECKS(s>=0,"surface name in surface_mol is not recognized"); }
981 ps=surfstring2ps(shapenm);
982 CHECKS(ps!=PSnone,"in surface_mol, panel shape name not recognized");
983 line2=strnword(line2,4);
984 if(line2) {
985 for(d=0;d<dim;d++) {
986 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&v1[d]);
987 CHECKS(itct==1,"incomplete vector in drift");
988 line2=strnword(line2,2); }
989 CHECKS(s>=0 && ps!=PSall && strcmp(nm1,"all"),"in surface_mol, use of coordinates requires that a specific panel be specified");
990 er=addsurfmol(sim,nmol,i,ms,v1,NULL,s,ps,nm1);
991 CHECKS(er!=1,"in surface_mol, unable to allocate temporary storage space");
992 CHECKS(er!=2,"in surface_mol, panel name not recognized");
993 CHECKS(er!=3,"not enough molecules permitted by max_mol"); }
994 else {
995 er=addsurfmol(sim,nmol,i,ms,NULL,NULL,s,ps,nm1);
996 CHECKS(er!=1,"in surface_mol, unable to allocate temporary storage space");
997 CHECKS(er!=2,"in surface_mol, no panels match the given description");
998 CHECKS(er!=3,"not enough molecules permitted by max_mol"); }
999 CHECKS(!line2,"unexpected text following surface_mol"); }
1000
1001 else if(!strcmp(word,"compartment_mol")) { // compartment_mol
1002 CHECKS(sim->mols,"need to enter species before compartment_mol");
1003 CHECKS(sim->cmptss,"compartments need to be defined before compartment_mol statement");
1004 CHECKS(sim->cmptss->ncmpt>0,"at least one compartment needs to be defined before compartment_mol");
1005 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&nmol);
1006 CHECKS(itct==1,"compartment_mol format: nmol species compartment");
1007 CHECKS(nmol>=0,"the number of molecules needs to be at least 0");
1008 CHECKS(line2=strnword(line2,2),"compartment_mol format: nmol species compartment");
1009 i=molstring2index1(sim,line2,&ms,NULL);
1010 CHECKS(i!=0,"wildcards and species groups are not permitted");
1011 CHECKS(i!=-1,"error reading molecule name");
1012 CHECKS(i!=-2,"mismatched or improper parentheses around molecule state");
1013 CHECKS(i!=-3,"cannot read molecule state value");
1014 CHECKS(i!=-4,"molecule name not recognized");
1015 CHECKS(i!=-5,"'all' species is not permitted");
1016 CHECKS(i!=-6,"cannot read wildcard logic");
1017 CHECKS(i!=-7,"out of memory");
1018 CHECKS(ms==MSsoln,"state needs to be solution");
1019 CHECKS(line2=strnword(line2,2),"compartment_mol format: nmol species compartment");
1020 itct=sscanf(line2,"%s",nm);
1021 CHECKS(itct==1,"compartment_mol format: nmol species compartment");
1022 c=stringfind(sim->cmptss->cnames,sim->cmptss->ncmpt,nm);
1023 CHECKS(c>=0,"compartment name is not recognized");
1024 er=addcompartmol(sim,nmol,i,sim->cmptss->cmptlist[c]);
1025 CHECKS(er!=2,"compartment volume is zero or nearly zero");
1026 CHECKS(er!=3,"not enough molecules permitted by max_mol");
1027 CHECKS(!strnword(line2,2),"unexpected text following compartment_mol"); }
1028
1029 else if(!strcmp(word,"molecule_lists")) { // molecule_lists
1030 while(line2) {
1031 itct=sscanf(line2,"%s",nm);
1032 CHECKS(itct==1,"unable to read molecule list name");
1033 er=addmollist(sim,nm,MLTsystem);
1034 CHECKS(er!=-1,"SMOLDYN BUG: out of memory");
1035 CHECKS(er!=-2,"molecule list name has already been used");
1036 CHECKS(er!=-3,"SMOLDYN BUG: illegal addmollist inputs");
1037 line2=strnword(line2,2); }
1038 CHECKS(!line2,"unexpected text following molecule_lists"); }
1039
1040 else if(!strcmp(word,"mol_list")) { // mol_list
1041 CHECKS(sim->mols && sim->mols->nlist>0,"need to enter molecule_lists before mol_list");
1042 er=molstring2index1(sim,line2,&ms,&index);
1043 CHECKS(er!=-1,"mol_list format: species[(state)] list_name");
1044 CHECKS(er!=-2,"mismatched or improper parentheses around molecule state");
1045 CHECKS(er!=-3,"cannot read molecule state value");
1046 CHECKS(er!=-4,"molecule name not recognized");
1047 CHECKS(er!=-6,"cannot read wildcard logic");
1048 CHECKS(er!=-7,"error allocating memory");
1049 CHECKS(ms<MSMAX || ms==MSall,"invalid state");
1050 CHECKS(line2=strnword(line2,2),"missing list name for mol_list");
1051 itct=sscanf(line2,"%s",nm);
1052 CHECKS(itct==1,"mol_list format: name[(state)] list_name");
1053 ll=stringfind(sim->mols->listname,sim->mols->nlist,nm);
1054 CHECKS(ll>=0,"molecule list name is not recognized");
1055 molsetlistlookup(sim,0,index,ms,ll);
1056 CHECKS(!strnword(line2,2),"unexpected text following mol_list"); }
1057
1058 else if(!strcmp(word,"mol_list_rule")) { // mol_list_rule
1059 CHECKS(sim->mols && sim->mols->nlist>0,"need to enter molecule_lists before mol_list");
1060 er=molstring2pattern(line2,&ms,pattern,0);
1061 CHECKS(er!=-1,"mol_list_pattern format: species[(state)] list_name");
1062 CHECKS(er!=-2,"mismatched or improper parentheses around molecule state");
1063 CHECKS(er!=-3,"cannot read molecule state value");
1064 CHECKS(ms<MSMAX || ms==MSall,"invalid state");
1065 CHECKS(line2=strnword(line2,2),"missing list name for mol_list");
1066 itct=sscanf(line2,"%s",nm);
1067 CHECKS(itct==1,"mol_list_pattern format: name[(state)] list_name");
1068 ll=stringfind(sim->mols->listname,sim->mols->nlist,nm);
1069 CHECKS(ll>=0,"molecule list name is not recognized");
1070 detailsi[0]=ll;
1071 er=RuleAddRule(sim,RTmollist,NULL,pattern,&ms,NULL,0,detailsi,NULL);
1072 CHECKS(!er,"out of memory in mol_list_rule");
1073 CHECKS(!strnword(line2,2),"unexpected text following mol_list_rule"); }
1074
1075 // graphics
1076
1077 else if(!strcmp(word,"graphics")) { // graphics
1078 itct=sscanf(line2,"%s",nm);
1079 CHECKS(itct==1,"missing graphics parameter");
1080 er=graphicsenablegraphics(sim,nm);
1081 CHECKS(er!=1,"out of memory");
1082 CHECKS(er!=2,"BUG: missing parameter");
1083 CHECKS(er!=3,"graphics method not recognized");
1084 CHECKS(!strnword(line2,2),"unexpected text following graphics"); }
1085
1086 else if(!strcmp(word,"graphic_iter")) { // graphic_iter
1087 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&i1);
1088 CHECKS(itct==1,"graphic_iter need to be a number");
1089 er=graphicssetiter(sim,i1);
1090 CHECKS(er!=1,"out of memory enabling graphics");
1091 CHECKS(er!=2,"BUG: missing parameter");
1092 CHECKS(er!=3,"graphic_iter needs to be >=1");
1093 CHECKS(!strnword(line2,2),"unexpected text following graphic_iter"); }
1094
1095 else if(!strcmp(word,"graphic_delay")) { // graphic_delay
1096 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&i1);
1097 CHECKS(itct==1,"graphic_delay need to be a number");
1098 er=graphicssetdelay(sim,i1);
1099 CHECKS(er!=1,"out of memory enabling graphics");
1100 CHECKS(er!=2,"BUG: missing parameter");
1101 CHECKS(er!=3,"graphic_delay needs to be >=0");
1102 CHECKS(!strnword(line2,2),"unexpected text following graphic_delay"); }
1103
1104 else if(!strcmp(word,"quit_at_end")) { // quit_at_end
1105 itct=sscanf(line2,"%s",nm);
1106 CHECKS(itct==1,"quit_at_end format: yes or no");
1107 if(strbegin(nm,"yes",0) || !strcmp(nm,"1")) sim->quitatend=1;
1108 else if(strbegin(nm,"no",0) || !strcmp(nm,"0")) sim->quitatend=0;
1109 else CHECKS(0,"cannot read quit_at_end value");
1110 CHECKS(!strnword(line2,2),"unexpected text following quit_at_end"); }
1111
1112 else if(!strcmp(word,"frame_thickness")) { // frame_thickness
1113 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
1114 CHECKS(itct==1,"frame_thickness needs to be a number");
1115 er=graphicssetframethickness(sim,flt1);
1116 CHECKS(er!=1,"out of memory enabling graphics");
1117 CHECKS(er!=2,"BUG: missing parameter");
1118 CHECKS(er!=3,"frame_thickness needs to be >=0");
1119 CHECKS(!strnword(line2,2),"unexpected text following frame_thickness"); }
1120
1121 else if(!strcmp(word,"frame_color") || !strcmp(word,"frame_colour")) { // frame_color
1122 er=graphicsreadcolor(&line2,v2);
1123 CHECKS(er!=3,"color values need to be between 0 and 1");
1124 CHECKS(er!=4,"color name not recognized");
1125 CHECKS(er!=6,"alpha values need to be between 0 and 1");
1126 CHECKS(er==0,"format is either 3 numbers or color name, and then optional alpha value");
1127 er=graphicssetframecolor(sim,v2);
1128 CHECKS(er!=1,"out of memory enabling graphics");
1129 CHECKS(er!=2,"BUG: missing parameter");
1130 CHECKS(er!=3,"frame_color values need to be between 0 and 1");
1131 CHECKS(!line2,"unexpected text following frame_color"); }
1132
1133 else if(!strcmp(word,"grid_thickness")) { // grid_thickness
1134 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
1135 CHECKS(itct==1,"grid_thickness needs to be a number");
1136 er=graphicssetgridthickness(sim,flt1);
1137 CHECKS(er!=1,"out of memory enabling graphics");
1138 CHECKS(er!=2,"BUG: missing parameter");
1139 CHECKS(er!=3,"grid_thickness needs to be >=0");
1140 CHECKS(!strnword(line2,2),"unexpected text following grid_thickness"); }
1141
1142 else if(!strcmp(word,"grid_color") || !strcmp(word,"grid_colour")) { // grid_color
1143 er=graphicsreadcolor(&line2,v2);
1144 CHECKS(er!=3,"color values need to be between 0 and 1");
1145 CHECKS(er!=4,"color name not recognized");
1146 CHECKS(er!=6,"alpha values need to be between 0 and 1");
1147 CHECKS(er==0,"format is either 3 numbers or color name, and then optional alpha value");
1148 er=graphicssetgridcolor(sim,v2);
1149 CHECKS(er!=1,"out of memory enabling graphics");
1150 CHECKS(er!=2,"BUG: missing parameter");
1151 CHECKS(er!=3,"grid_color values need to be between 0 and 1");
1152 CHECKS(!line2,"unexpected text following grid_color"); }
1153
1154 else if(!strcmp(word,"background_color") || !strcmp(word,"background_colour")) { // background_color
1155 er=graphicsreadcolor(&line2,v2);
1156 CHECKS(er!=3,"color values need to be between 0 and 1");
1157 CHECKS(er!=4,"color name not recognized");
1158 CHECKS(er!=6,"alpha values need to be between 0 and 1");
1159 CHECKS(er==0,"format is either 3 numbers or color name, and then optional alpha value");
1160 er=graphicssetbackcolor(sim,v2);
1161 CHECKS(er!=1,"out of memory enabling graphics");
1162 CHECKS(er!=2,"BUG: missing parameter");
1163 CHECKS(er!=3,"background_color values need to be between 0 and 1");
1164 CHECKS(!line2,"unexpected text following background_color"); }
1165
1166 else if(!strcmp(word,"text_color") || !strcmp(word,"text_colour")) { // text_color
1167 er=graphicsreadcolor(&line2,v2);
1168 CHECKS(er!=3,"color values need to be between 0 and 1");
1169 CHECKS(er!=4,"color name not recognized");
1170 CHECKS(er!=6,"alpha values need to be between 0 and 1");
1171 CHECKS(er==0,"format is either 3 numbers or color name, and then optional alpha value");
1172 er=graphicssetbackcolor(sim,v2);
1173 CHECKS(er!=1,"out of memory enabling graphics");
1174 CHECKS(er!=2,"BUG: missing parameter");
1175 CHECKS(er!=3,"text_color values need to be between 0 and 1");
1176 CHECKS(!line2,"unexpected text following text_color"); }
1177
1178 else if(!strcmp(word,"text_display")) { // text_display
1179 while(line2) {
1180 itct=sscanf(line2,"%s",nm);
1181 CHECKS(itct==1,"error reading text_display item");
1182 er=graphicssettextitem(sim,nm);
1183 CHECKS(er!=1,"out of memory");
1184 CHECKS(er!=2,"unrecognized text display item (check that species have been defined)");
1185 line2=strnword(line2,2); }
1186 CHECKS(!line2,"unexpected text following text_display"); }
1187
1188 else if(!strcmp(word,"light")) { // light
1189 itct=sscanf(line2,"%s",nm);
1190 CHECKS(itct==1,"light format: light_number parameter values");
1191 if(!strcmp(nm,"global") || !strcmp(nm,"room")) lt=-1;
1192 else {
1193 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,<);
1194 CHECKS(lt>=0 && lt<MAXLIGHTS,"light number out of bounds"); }
1195 line2=strnword(line2,2);
1196 CHECKS(line2,"light format: light_number parameter values");
1197 itct=sscanf(line2,"%s",nm);
1198 CHECKS(itct==1,"unable to read lighting parameter");
1199 ltparam=graphicsstring2lp(nm);
1200 CHECKS(ltparam!=LPnone,"unrecognized light parameter");
1201 line2=strnword(line2,2);
1202 if(ltparam==LPon || ltparam==LPoff || ltparam==LPauto) v2[0]=0;
1203 else {
1204 CHECKS(line2,"light format: light_number parameter values");
1205 itct=strmathsscanf(line2,"%mlg %mlg %mlg",varnames,varvalues,nvar,&v2[0],&v2[1],&v2[2]);
1206 CHECKS(itct==3,"light is missing one or more values");
1207 v2[3]=1;
1208 line2=strnword(line2,4);
1209 if(line2) {
1210 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&v2[3]);
1211 CHECKS(itct==1,"failed to read alpha light value");
1212 line2=strnword(line2,2); }
1213 if(ltparam!=LPposition) {
1214 for(i1=0;i1<4;i1++)
1215 CHECKS(v2[i1]>=0 && v2[i1]<=1,"light color values need to be between 0 and 1"); }}
1216 er=graphicssetlight(sim,NULL,lt,ltparam,v2);
1217 CHECKS(er!=1,"out of memory");
1218 CHECKS(!er,"BUG: error in light statement");
1219 CHECKS(!line2,"unexpected text following light"); }
1220
1221 else if(!strcmp(word,"display_size")) { // display_size
1222 CHECKS(sim->mols,"need to enter species before display_size");
1223 er=molstring2index1(sim,line2,&ms,&index);
1224 CHECKS(er!=-1,"display_size format: species[(state)] value");
1225 CHECKS(er!=-2,"mismatched or improper parentheses around molecule state");
1226 CHECKS(er!=-3,"cannot read molecule state value");
1227 CHECKS(er!=-4,"molecule name not recognized");
1228 CHECKS(er!=-6,"cannot read wildcard logic");
1229 CHECKS(er!=-7,"error allocating memory");
1230 CHECKS(ms<MSMAX || ms==MSall,"invalid state");
1231 CHECKS(line2=strnword(line2,2),"missing data for display_size");
1232 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
1233 CHECKS(itct==1,"display_size format: name[(state)] size");
1234 CHECKS(flt1>=0,"display_size value needs to be >=0");
1235 molsetdisplaysize(sim,0,index,ms,flt1);
1236 CHECKS(!strnword(line2,2),"unexpected text following display_size"); }
1237
1238 else if(!strcmp(word,"display_size_rule")) { // display_size_rule
1239 er=molstring2pattern(line2,&ms,pattern,0);
1240 CHECKS(er!=-1,"display_size_rule format: species[(state)] value");
1241 CHECKS(er!=-2,"mismatched or improper parentheses around molecule state");
1242 CHECKS(er!=-3,"cannot read molecule state value");
1243 CHECKS(ms<MSMAX || ms==MSall,"invalid state");
1244 CHECKS(line2=strnword(line2,2),"missing data for display_size");
1245 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
1246 CHECKS(itct==1,"display_size_rule format: name[(state)] size");
1247 CHECKS(flt1>=0,"display_size_rule value needs to be >=0");
1248 er=RuleAddRule(sim,RTdispsize,NULL,pattern,&ms,NULL,flt1,NULL,NULL);
1249 CHECKS(!er,"out of memory in display_size_rule");
1250 CHECKS(!strnword(line2,2),"unexpected text following display_size_rule"); }
1251
1252 else if(!strcmp(word,"color") || !strcmp(word,"colour")) { // color
1253 CHECKS(sim->mols,"need to enter species before color");
1254 er=molstring2index1(sim,line2,&ms,&index);
1255 CHECKS(er!=-1,"color format: species[(state)] color");
1256 CHECKS(er!=-2,"mismatched or improper parentheses around molecule state");
1257 CHECKS(er!=-3,"cannot read molecule state value");
1258 CHECKS(er!=-4,"molecule name not recognized");
1259 CHECKS(er!=-6,"cannot read wildcard logic");
1260 CHECKS(er!=-7,"error allocating memory");
1261 CHECKS(ms<MSMAX || ms==MSall,"invalid state");
1262 line2=strnword(line2,2);
1263 er=graphicsreadcolor(&line2,v2);
1264 CHECKS(er!=3,"color values need to be between 0 and 1");
1265 CHECKS(er!=4,"color name not recognized");
1266 CHECKS(er!=6,"alpha values need to be between 0 and 1");
1267 CHECKS(er==0,"format is either 3 numbers or color name, and then optional alpha value");
1268 molsetcolor(sim,0,index,ms,v2);
1269 CHECKS(!line2,"unexpected text following color"); }
1270
1271 else if(!strcmp(word,"color_rule") || !strcmp(word,"colour_rule")) { // color_rule
1272 er=molstring2pattern(line2,&ms,pattern,0);
1273 CHECKS(er!=-1,"color format: species[(state)] color");
1274 CHECKS(er!=-2,"mismatched or improper parentheses around molecule state");
1275 CHECKS(er!=-3,"cannot read molecule state value");
1276 CHECKS(ms<MSMAX || ms==MSall,"invalid state");
1277 line2=strnword(line2,2);
1278 er=graphicsreadcolor(&line2,v2);
1279 CHECKS(er!=3,"color values need to be between 0 and 1");
1280 CHECKS(er!=4,"color name not recognized");
1281 CHECKS(er!=6,"alpha values need to be between 0 and 1");
1282 CHECKS(er==0,"format is either 3 numbers or color name, and then optional alpha value");
1283 er=RuleAddRule(sim,RTcolor,NULL,pattern,&ms,NULL,0,NULL,v2);
1284 CHECKS(!er,"out of memory in color_rule");
1285 CHECKS(!line2,"unexpected text following color_rule"); }
1286
1287 else if(!strcmp(word,"tiff_iter")) { // tiff_iter
1288 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&i1);
1289 CHECKS(itct==1,"tiff_iter needs to be a number");
1290 er=graphicssettiffiter(sim,i1);
1291 CHECKS(er!=1,"out of memory enabling graphics");
1292 CHECKS(er!=2,"BUG: missing parameter");
1293 CHECKS(er!=3,"tiff_iter needs to be >=1");
1294 CHECKS(!strnword(line2,2),"unexpected text following tiff_iter"); }
1295
1296 else if(!strcmp(word,"tiff_name")) { // tiff_name
1297 itct=sscanf(line2,"%s",nm);
1298 CHECKS(itct==1,"tiff_name needs to be a string");
1299 strcpy(nm1,sim->filepath);
1300 CHECKS(strlen(nm)<STRCHAR-strlen(nm1),"tiff name is too long");
1301 strcat(nm1,nm);
1302 gl2SetOptionStr("TiffName",nm1);
1303 CHECKS(!strnword(line2,2),"unexpected text following tiff_name"); }
1304
1305 else if(!strcmp(word,"tiff_min")) { // tiff_min
1306 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&i1);
1307 CHECKS(itct==1,"tiff_min needs to be a number");
1308 CHECKS(i1>=0,"tiff_min has to be at least 0");
1309 gl2SetOptionInt("TiffNumber",i1);
1310 CHECKS(!strnword(line2,2),"unexpected text following tiff_min"); }
1311
1312 else if(!strcmp(word,"tiff_max")) { // tiff_max
1313 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&i1);
1314 CHECKS(itct==1,"tiff_max needs to be a number");
1315 CHECKS(i1>=0,"tiff_max has to be at least 0");
1316 gl2SetOptionInt("TiffNumMax",i1);
1317 CHECKS(!strnword(line2,2),"unexpected text following tiff_max"); }
1318
1319 // about runtime commands
1320
1321 else if(!strcmp(word,"output_root")) { // output_root
1322 er=scmdsetfroot(sim->cmds,line2);
1323 CHECKS(er!=-1,"SMOLDYN BUG: scmdsetfroot"); }
1324
1325 else if(!strcmp(word,"output_files") || !strcmp(word,"output_file")) { // output_files
1326 er=scmdsetfnames(sim->cmds,line2,0);
1327 CHECKS(er!=1,"out of memory in output_files");
1328 CHECKS(er!=2,"error reading file name");
1329 CHECKS(er!=4,"BUG: variable cmds became NULL"); }
1330
1331 else if(!strcmp(word,"output_data")) { // output_data
1332 er=scmdsetdnames(sim->cmds,line2);
1333 CHECKS(er!=1,"out of memory in output_data");
1334 CHECKS(er!=2,"error reading data variable name");
1335 CHECKS(er!=4,"BUG: variable cmds became NULL"); }
1336
1337 else if(!strcmp(word,"output_precision")) { // output_precision
1338 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&i1);
1339 CHECKS(itct==1,"format for output_precision: value");
1340 scmdsetprecision(sim->cmds,i1);
1341 CHECKS(!strnword(line2,2),"unexpected text following output_precision"); }
1342
1343 else if(!strcmp(word,"output_format")) { // output_format
1344 itct=sscanf(line2,"%s",nm);
1345 CHECKS(itct==1,"format for output_format: csv or ssv");
1346 er=scmdsetoutputformat(sim->cmds,nm);
1347 CHECKS(!er,"output_format value not recognized");
1348 CHECKS(!strnword(line2,2),"unexpected text following output_format"); }
1349
1350 else if(!strcmp(word,"append_files")) { // append_files
1351 er=scmdsetfnames(sim->cmds,line2,1);
1352 CHECKS(er!=1,"out of memory");
1353 CHECKS(er!=2,"error reading file name");
1354 CHECKS(er!=4,"BUG: variable cmds became NULL"); }
1355
1356 else if(!strcmp(word,"output_file_number")) { // output_file_number
1357 itct=strmathsscanf(line2,"%s %mi",varnames,varvalues,nvar,nm,&i1);
1358 CHECKS(itct==2,"format for output_file_number: filename number");
1359 CHECKS(i1>=0,"output_file_number needs to be >= 0");
1360 er=scmdsetfsuffix(sim->cmds,nm,i1);
1361 CHECKS(!er,"error setting output_file_number");
1362 CHECKS(!strnword(line2,3),"unexpected text following output_file_number"); }
1363
1364 else if(!strcmp(word,"cmd")) { // cmd
1365 er=scmdstr2cmd(sim->cmds,line2,varnames,varvalues,nvar);
1366 CHECKS(er!=1,"out of memory in cmd");
1367 CHECKS(er!=2,"BUG: no command superstructure for cmd");
1368 CHECKS(er!=3,"cmd format: type [on off dt] string");
1369 CHECKS(er!=6,"command timing type character not recognized"); }
1370
1371 else if(!strcmp(word,"max_cmd")) { // max_cmd
1372 CHECKS(0,"max_cmd is an obsolete statement. It simply needs to be removed."); }
1373
1374 // surfaces
1375
1376 else if(!strcmp(word,"max_surface")) { // max_surface
1377 CHECKS(dim>0,"need to enter dim before max_surface");
1378 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&i1);
1379 CHECKS(itct==1,"max_surface needs to be a number");
1380 CHECKS(i1>0,"max_surface must be greater than 0");
1381 CHECKS(surfenablesurfaces(sim,i1)==0,"failed to enable surfaces");
1382 CHECKS(!strnword(line2,2),"unexpected text following max_surface"); }
1383
1384 else if(!strcmp(word,"new_surface")) { // new_surface
1385 CHECKS(dim>0,"need to enter dim before new_surface");
1386 srf=surfreadstring(sim,pfp,NULL,"name",line2);
1387 CHECK(srf!=NULL); }
1388
1389 else if(!strcmp(word,"surface")) { // surface
1390 CHECKS(dim>0,"need to enter dim before surface");
1391 CHECKS(sim->srfss,"individual surfaces need to be defined before using surface");
1392 itct=sscanf(line2,"%s %s",nm,nm1);
1393 CHECKS(itct==2,"surface format: surface_name statement_name statement_text");
1394 line2=strnword(line2,3);
1395 CHECKS(line2,"surface format: surface_name statement_name statement_text");
1396 s=stringfind(sim->srfss->snames,sim->srfss->nsrf,nm);
1397 CHECKS(s>=0,"surface is unrecognized");
1398 srf=sim->srfss->srflist[s];
1399 srf=surfreadstring(sim,pfp,srf,nm1,line2);
1400 CHECK(srf!=NULL); }
1401
1402 // compartments
1403
1404 else if(!strcmp(word,"max_compartment")) { // max_compartment
1405 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&i1);
1406 CHECKS(itct==1,"max_compartment needs to be a number");
1407 CHECKS(i1>=0,"max_compartment must be at least 0");
1408 CHECKS(!compartenablecomparts(sim,i1),"failed to enable compartments");
1409 CHECKS(!strnword(line2,2),"unexpected text following max_compartment"); }
1410
1411 else if(!strcmp(word,"new_compartment")) { // new_compartment
1412 cmpt=compartreadstring(sim,pfp,NULL,"name",line2);
1413 CHECK(cmpt!=NULL); }
1414
1415 else if(!strcmp(word,"compartment")) { // compartment
1416 CHECKS(sim->cmptss,"individual compartments need to be defined before using compartment");
1417 itct=sscanf(line2,"%s %s",nm,nm1);
1418 CHECKS(itct==2,"compartment format: compart_name statement_name statement_text");
1419 line2=strnword(line2,3);
1420 CHECKS(line2,"compartment format: compart_name statement_name statement_text");
1421 c=stringfind(sim->cmptss->cnames,sim->cmptss->ncmpt,nm);
1422 CHECKS(c>=0,"compartment is unrecognized");
1423 cmpt=sim->cmptss->cmptlist[c];
1424 cmpt=compartreadstring(sim,pfp,cmpt,nm1,line2);
1425 CHECK(cmpt!=NULL); }
1426
1427 // ports
1428
1429 else if(!strcmp(word,"max_port")) { // max_port
1430 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&i1);
1431 CHECKS(itct==1,"max_port needs to be a number");
1432 CHECKS(i1>=0,"max_port must be at least 0");
1433 CHECKS(portenableports(sim,i1),"failed to enable ports");
1434 CHECKS(!strnword(line2,2),"unexpected text following max_port"); }
1435
1436 else if(!strcmp(word,"new_port")) { // new_port
1437 port=portreadstring(sim,pfp,NULL,"name",line2);
1438 CHECK(port!=NULL); }
1439
1440 else if(!strcmp(word,"port")) { // port
1441 CHECKS(sim->portss,"individual ports need to be defined before using port");
1442 itct=sscanf(line2,"%s %s",nm,nm1);
1443 CHECKS(itct==2,"port format: port_name statement_name statement_text");
1444 line2=strnword(line2,3);
1445 CHECKS(line2,"port format: port_name statement_name statement_text");
1446 prt=stringfind(sim->portss->portnames,sim->portss->nport,nm);
1447 CHECKS(prt>=0,"port is unrecognized");
1448 port=sim->portss->portlist[prt];
1449 port=portreadstring(sim,pfp,port,nm1,line2);
1450 CHECK(port!=NULL); }
1451
1452 // filaments
1453
1454 else if(!strcmp(word,"max_filament")) { // max_filament
1455 CHECKS(0,"max_filament has been deprecated"); }
1456
1457 /*
1458 else if(!strcmp(word,"new_filament")) { // new_filament
1459 fil=filreadstring(sim,pfp,NULL,"name",line2);
1460 CHECK(fil!=NULL); }
1461
1462 else if(!strcmp(word,"filament")) { // filament
1463 CHECKS(sim->filss,"individual filaments need to be defined before using filament");
1464 itct=sscanf(line2,"%s %s",nm,nm1);
1465 CHECKS(itct==2,"filament format: filament_name statement_name statement_text");
1466 line2=strnword(line2,3);
1467 CHECKS(line2,"filament format: filament_name statement_name statement_text");
1468 f=stringfind(sim->filss->fnames,sim->filss->nfil,nm);
1469 CHECKS(f>=0,"filament is unrecognized");
1470 fil=sim->filss->fillist[f];
1471 fil=filreadstring(sim,pfp,fil,nm1,line2);
1472 CHECK(fil!=NULL); }
1473 */
1474
1475 else if(!strcmp(word,"random_filament")) { // random_filament
1476 CHECKS(sim->filss,"need to enter a filament type before random_filament");
1477 filss=sim->filss;
1478 itct=sscanf(line2,"%s %s",nm,nm1);
1479 CHECKS(itct==2,"random_filament format: name type segments [x y z]");
1480 ft=stringfind(filss->ftnames,filss->ntype,nm1);
1481 CHECKS(ft>=0,"filament type is unknown");
1482 filtype=filss->filtypes[ft];
1483 fil=filAddFilament(filtype,NULL,nm);
1484 CHECKS(fil,"unable to add filament to simulation");
1485 line2=strnword(line2,3);
1486
1487 CHECKS(line2,"random_filament format: name type segments [x y z]");
1488 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&i1);
1489 CHECKS(itct==1,"random_segments format: number [x y z] [thickness]");
1490 CHECKS(i1>0,"number needs to be >0");
1491 line2=strnword(line2,2);
1492 if(fil->nbs==0) {
1493 CHECKS(line2,"missing x y z position information");
1494 itct=sscanf(line2,"%s %s %s",str1,str2,str3);
1495 CHECKS(itct==3,"random_segments format: number [x y z] [thickness]");
1496 line2=strnword(line2,4); }
1497 else {
1498 sprintf(str1,"%i",0);
1499 sprintf(str2,"%i",0);
1500 sprintf(str3,"%i",0); }
1501 thick=1;
1502 if(line2) {
1503 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&thick);
1504 CHECKS(itct==1,"random_segments format: number [x y z] [thickness]");
1505 CHECKS(thick>0,"thickness needs to be >0");
1506 line2=strnword(line2,2); }
1507 if(filtype->isbead)
1508 er=filAddRandomBeads(fil,i1,str1,str2,str3);
1509 else
1510 er=filAddRandomSegments(fil,i1,str1,str2,str3,thick);
1511 CHECKS(er!=2,"random_segments positions need to be 'u' or value");
1512 CHECKS(er==0,"BUG: error in filAddRandomsegments");
1513 CHECKS(!line2,"unexpected text following random_segments"); }
1514
1515
1516 // reactions
1517
1518 else if(!strcmp(word,"reaction") || !strcmp(word,"reaction_rule") || !strcmp(word,"reaction_cmpt") || !strcmp(word,"reaction_surface")) { // reaction, reaction_rule, reaction_cmpt, reaction_surface
1519 CHECKS(sim->mols,"need to enter species before reaction");
1520 CHECKS(sim->dim>0,"need to enter system dimensionality before reaction");
1521 er=rxnparsereaction(sim,word,line2,errstr);
1522 CHECKS(!er,"%s",errstr); }
1523
1524 #ifdef OPTION_VCELL
1525 else if(!strcmp(word,"reaction_rate")) { // reaction_rate
1526 if(line2){
1527 stringstream ss(line2);
1528 ss >> rname;
1529 string rawStr, expStr;
1530 getline(ss, rawStr);
1531 size_t found = rawStr.find(";");
1532
1533 if(found!=string::npos) {
1534 expStr = rawStr.substr(0, found);
1535 ValueProvider* valueProvider = sim->valueProviderFactory->createValueProvider(expStr);
1536 double constRate = 1;
1537 try {
1538 constRate = valueProvider->getConstantValue();
1539 } catch (...) {
1540 rxn->rateValueProvider = valueProvider;
1541 }
1542
1543 //if rate is a exp, a "fake" set value has to be set to allow RxnSetRate to do proper job, especially when there
1544 //is reversible reactions, the rxn->rparamt, and rxn->bindrad2 have to be set.
1545 er = RxnSetValue(sim, "rate", rxn, constRate);
1546 line2=NULL; }
1547 else {
1548 itct=sscanf(line2,"%s %lg",rname,&flt1);
1549 CHECKS(itct==2,"reaction_rate format: rname rate");
1550 r=readrxnname(sim,rname,&order,&rxn,NULL,1);
1551 CHECKS(r>=0,"unrecognized reaction name");
1552 er=RxnSetValue(sim,"rate",rxn,flt1);
1553 CHECKS(er!=4,"reaction rate value must be non-negative");
1554 CHECKS(!strnword(line2,3),"unexpected text following reaction"); }}}
1555 #else
1556
1557 else if(!strcmp(word,"reaction_rate")) { // reaction_rate
1558 itct=strmathsscanf(line2,"%s %mlg",varnames,varvalues,nvar,rname,&flt1);
1559 CHECKS(itct==2,"reaction_rate format: rname rate");
1560 CHECKS(flt1>=0,"reaction rate value must be non-negative");
1561 rxn=NULL;
1562 r=readrxnname(sim,rname,NULL,&rxn,NULL,1);
1563 if(r>=0)
1564 RxnSetValue(sim,"rate",rxn,flt1);
1565 r=readrxnname(sim,rname,NULL,&rxn,&vlist,2);
1566 CHECKS(r!=-2,"out of memory reading reaction name");
1567 if(r>=0) {
1568 for(i=0;i<vlist->n;i++)
1569 RxnSetValue(sim,"rate",(rxnptr) vlist->xs[i],flt1);
1570 ListFreeV(vlist); }
1571 r=readrxnname(sim,rname,NULL,&rxn,NULL,3);
1572 CHECKS(r!=-2,"out of memory reading reaction name");
1573 if(r>=0)
1574 RxnSetValue(NULL,"rate",rxn,flt1);
1575 CHECKS(rxn,"unrecognized reaction name");
1576 CHECKS(!strnword(line2,3),"unexpected text following reaction_rate"); }
1577 #endif
1578
1579 else if(!strcmp(word,"reaction_multiplicity")) { // reaction_multiplicity
1580 itct=strmathsscanf(line2,"%s %mi",varnames,varvalues,nvar,rname,&i1);
1581 CHECKS(itct==2,"reaction_multiplicity format: rname multiplicity");
1582 CHECKS(i1>=0,"reaction multiplicity value must be non-negative");
1583 rxn=NULL;
1584 r=readrxnname(sim,rname,NULL,&rxn,NULL,1);
1585 if(r>=0)
1586 RxnSetValue(sim,"multiplicity",rxn,(double)i1);
1587 r=readrxnname(sim,rname,NULL,&rxn,&vlist,2);
1588 CHECKS(r!=-2,"out of memory reading reaction name");
1589 if(r>=0) {
1590 for(i=0;i<vlist->n;i++)
1591 RxnSetValue(sim,"multiplicity",(rxnptr) vlist->xs[i],(double)i1);
1592 ListFreeV(vlist); }
1593 r=readrxnname(sim,rname,NULL,&rxn,NULL,3);
1594 CHECKS(r!=-2,"out of memory reading reaction name");
1595 if(r>=0)
1596 RxnSetValue(NULL,"multiplicity",rxn,(double)i1);
1597 CHECKS(rxn,"unrecognized reaction name");
1598 CHECKS(!strnword(line2,3),"unexpected text following reaction_multiplicity"); }
1599
1600 else if(!strcmp(word,"confspread_radius")) { // confspread_radius
1601 itct=strmathsscanf(line2,"%s %mlg",varnames,varvalues,nvar,rname,&flt1);
1602 CHECKS(itct==2,"confspread_radius format: rname radius");
1603 CHECKS(flt1>=0,"confspread radius value must be non-negative");
1604 rxn=NULL;
1605 r=readrxnname(sim,rname,NULL,&rxn,NULL,1);
1606 if(r>=0)
1607 RxnSetValue(sim,"confspreadrad",rxn,flt1);
1608 r=readrxnname(sim,rname,NULL,&rxn,&vlist,2);
1609 CHECKS(r!=-2,"out of memory reading reaction name");
1610 if(r>=0) {
1611 for(i=0;i<vlist->n;i++)
1612 RxnSetValue(sim,"confspreadrad",(rxnptr) vlist->xs[i],flt1);
1613 ListFreeV(vlist); }
1614 r=readrxnname(sim,rname,NULL,&rxn,NULL,3);
1615 CHECKS(r!=-2,"out of memory reading reaction name");
1616 if(r>=0)
1617 RxnSetValue(NULL,"confspreadrad",rxn,flt1);
1618 CHECKS(rxn,"unrecognized reaction name");
1619 CHECKS(!strnword(line2,3),"unexpected text following confspread_radius"); }
1620
1621 else if(!strcmp(word,"binding_radius")) { // binding_radius
1622 itct=strmathsscanf(line2,"%s %mlg",varnames,varvalues,nvar,rname,&flt1);
1623 CHECKS(itct==2,"binding_radius format: rname radius");
1624 CHECKS(flt1>=0,"binding radius value must be non-negative");
1625 rxn=NULL;
1626 r=readrxnname(sim,rname,NULL,&rxn,NULL,1);
1627 if(r>=0)
1628 RxnSetValue(sim,"bindrad",rxn,flt1);
1629 r=readrxnname(sim,rname,NULL,&rxn,&vlist,2);
1630 CHECKS(r!=-2,"out of memory reading reaction name");
1631 if(r>=0) {
1632 for(i=0;i<vlist->n;i++)
1633 RxnSetValue(sim,"bindrad",(rxnptr) vlist->xs[i],flt1);
1634 ListFreeV(vlist); }
1635 r=readrxnname(sim,rname,NULL,&rxn,NULL,3);
1636 CHECKS(r!=-2,"out of memory reading reaction name");
1637 if(r>=0)
1638 RxnSetValue(NULL,"bindrad",rxn,flt1);
1639 CHECKS(rxn,"unrecognized reaction name");
1640 CHECKS(!strnword(line2,3),"unexpected text following binding_radius"); }
1641
1642 else if(!strcmp(word,"reaction_probability")) { // reaction_probability
1643 itct=strmathsscanf(line2,"%s %mlg",varnames,varvalues,nvar,rname,&flt1);
1644 CHECKS(itct==2,"reaction_probability format: rname value");
1645 CHECKS(flt1>=0 && flt1<=1,"probability value must be between 0 and 1");
1646 rxn=NULL;
1647 r=readrxnname(sim,rname,NULL,&rxn,NULL,1);
1648 if(r>=0)
1649 RxnSetValue(sim,"prob",rxn,flt1);
1650 r=readrxnname(sim,rname,NULL,&rxn,&vlist,2);
1651 CHECKS(r!=-2,"out of memory reading reaction name");
1652 if(r>=0) {
1653 for(i=0;i<vlist->n;i++)
1654 RxnSetValue(sim,"prob",(rxnptr) vlist->xs[i],flt1);
1655 ListFreeV(vlist); }
1656 r=readrxnname(sim,rname,NULL,&rxn,NULL,3);
1657 CHECKS(r!=-2,"out of memory reading reaction name");
1658 if(r>=0)
1659 RxnSetValue(NULL,"prob",rxn,flt1);
1660 CHECKS(rxn,"unrecognized reaction name");
1661 CHECKS(!strnword(line2,3),"unexpected text following reaction_probability"); }
1662
1663 else if(!strcmp(word,"reaction_chi")) { // reaction_chi
1664 itct=strmathsscanf(line2,"%s %mlg",varnames,varvalues,nvar,rname,&flt1);
1665 CHECKS(itct==2,"reaction_chi format: rname value");
1666 CHECKS(flt1!=0 && flt1<1,"reaction chi value must be between 0 and 1 (or -1 to disable)");
1667 rxn=NULL;
1668 r=readrxnname(sim,rname,&order,&rxn,NULL,1);
1669 if(r>=0) {
1670 CHECKS(order==2,"reaction chi value is only allowed for order 2 reactions");
1671 RxnSetValue(sim,"chi",rxn,flt1); }
1672 r=readrxnname(sim,rname,&order,&rxn,&vlist,2);
1673 CHECKS(r!=-2,"out of memory reading reaction name");
1674 if(r>=0) {
1675 CHECKS(order==2,"reaction chi value is only allowed for order 2 reactions");
1676 for(i=0;i<vlist->n;i++)
1677 RxnSetValue(sim,"chi",(rxnptr) vlist->xs[i],flt1);
1678 ListFreeV(vlist); }
1679 r=readrxnname(sim,rname,&order,&rxn,NULL,3);
1680 CHECKS(r!=-2,"out of memory reading reaction name");
1681 if(r>=0) {
1682 CHECKS(order==2,"reaction chi value is only allowed for order 2 reactions");
1683 RxnSetValue(NULL,"chi",rxn,flt1); }
1684 CHECKS(rxn,"unrecognized reaction name");
1685 CHECKS(!strnword(line2,3),"unexpected text following reaction_chi"); }
1686
1687 else if(!strcmp(word,"reaction_production")) { // reaction_production
1688 itct=strmathsscanf(line2,"%s %mlg",varnames,varvalues,nvar,rname,&flt1);
1689 CHECKS(itct==2,"reaction_production format: rname value");
1690 CHECKS(flt1>=0 && flt1<=1,"production value must be between 0 and 1");
1691 rxn=NULL;
1692 r=readrxnname(sim,rname,&order,&rxn,NULL,1);
1693 if(r>=0) {
1694 CHECKS(order==0,"production is only allowed for order 0 reactions");
1695 RxnSetValue(sim,"prob",rxn,flt1); }
1696 r=readrxnname(sim,rname,&order,&rxn,&vlist,2);
1697 CHECKS(r!=-2,"out of memory reading reaction name");
1698 if(r>=0) {
1699 CHECKS(order==0,"production is only allowed for order 0 reactions");
1700 for(i=0;i<vlist->n;i++)
1701 RxnSetValue(sim,"prob",(rxnptr) vlist->xs[i],flt1);
1702 ListFreeV(vlist); }
1703 r=readrxnname(sim,rname,&order,&rxn,NULL,3);
1704 CHECKS(r!=-2,"out of memory reading reaction name");
1705 if(r>=0) {
1706 CHECKS(order==0,"production is only allowed for order 0 reactions");
1707 RxnSetValue(NULL,"prob",rxn,flt1); }
1708 CHECKS(rxn,"unrecognized reaction name");
1709 CHECKS(!strnword(line2,3),"unexpected text following reaction_production"); }
1710
1711 else if(!strcmp(word,"reaction_permit")) { // reaction_permit
1712 itct=sscanf(line2,"%s",rname);
1713 CHECKS(itct==1,"missing reaction name");
1714 r=readrxnname(sim,rname,&order,&rxn,NULL,1);
1715 CHECKS(r>=0,"unrecognized reaction name");
1716 CHECKS(order>0,"reaction_permit is not allowed for order 0 reactions");
1717 for(ord=0;ord<order;ord++) {
1718 CHECKS(line2=strnword(line2,2),"missing state term in reaction_permit");
1719 itct=sscanf(line2,"%s",nm1);
1720 rctstate[ord]=molstring2ms(nm1); }
1721 RxnSetPermit(sim,rxn,order,rctstate,1);
1722 CHECKS(!strnword(line2,3),"unexpected text following reaction_permit"); }
1723
1724 else if(!strcmp(word,"reaction_forbid")) { // reaction_forbid
1725 itct=sscanf(line2,"%s",rname);
1726 CHECKS(itct==1,"missing reaction name");
1727 r=readrxnname(sim,rname,&order,&rxn,&vlist,1);
1728 CHECKS(r>=0,"unrecognized reaction name");
1729 CHECKS(order>0,"reaction_forbid is not allowed for order 0 reactions");
1730 for(ord=0;ord<order;ord++) {
1731 CHECKS(line2=strnword(line2,2),"missing state term in reaction_forbid");
1732 itct=sscanf(line2,"%s",nm1);
1733 rctstate[ord]=molstring2ms(nm1); }
1734 RxnSetPermit(sim,rxn,order,rctstate,0);
1735 CHECKS(!strnword(line2,3),"unexpected text following reaction_forbid"); }
1736
1737 else if(!strcmp(word,"reaction_serialnum")) { // reaction_serialnum
1738 itct=sscanf(line2,"%s",rname);
1739 CHECKS(itct==1,"reaction_serialnum format: rname rules_list");
1740 rxn=NULL;
1741 readrxnname(sim,rname,&order,&rxn,NULL,1);
1742 if(!rxn) readrxnname(sim,rname,&order,&rxn,NULL,2);
1743 if(!rxn) readrxnname(sim,rname,&order,&rxn,NULL,3);
1744 CHECKS(rxn,"unrecognized reaction name");
1745 for(prd=0;prd<rxn->nprod;prd++) {
1746 sernolist[prd]=0;
1747 line2=strnword(line2,2);
1748 CHECKS(line2,"not enough product codes entered");
1749 itct=sscanf(line2,"%s",pattern);
1750 CHECKS(itct==1,"error reading a product code");
1751 if(!strcmp(pattern,"+")) prd--;
1752 else {
1753 pserno=rxnstring2sernocode(pattern,prd);
1754 CHECKS(pserno!=0,"error reading a product code");
1755 sernolist[prd]=pserno; }}
1756 rxn=NULL;
1757 r=readrxnname(sim,rname,NULL,&rxn,NULL,1);
1758 if(r>=0) {
1759 er=RxnSetPrdSerno(rxn,sernolist);
1760 CHECKS(!er,"out of memory allocating product serial number list"); }
1761 r=readrxnname(sim,rname,NULL,&rxn,&vlist,2);
1762 CHECKS(r!=-2,"out of memory reading reaction name");
1763 if(r>=0) {
1764 for(i=0;i<vlist->n;i++) {
1765 er=RxnSetPrdSerno((rxnptr) vlist->xs[i],sernolist);
1766 CHECKS(!er,"out of memory allocating product serial number list"); }
1767 ListFreeV(vlist); }
1768 r=readrxnname(sim,rname,NULL,&rxn,NULL,3);
1769 CHECKS(r!=-2,"out of memory reading reaction name");
1770 if(r>=0) {
1771 er=RxnSetPrdSerno(rxn,sernolist);
1772 CHECKS(!er,"out of memory allocating product serial number list"); }
1773 CHECKS(!strnword(line2,2),"unexpected text following reaction_serialnum"); }
1774
1775 else if(!strcmp(word,"reaction_intersurface")) { // reaction_intersurface
1776 itct=sscanf(line2,"%s",rname);
1777 CHECKS(itct==1,"reaction_intersurface format: rname rules_list");
1778 rxn=NULL;
1779 readrxnname(sim,rname,&order,&rxn,NULL,1);
1780 if(!rxn) readrxnname(sim,rname,&order,&rxn,NULL,2);
1781 if(!rxn) readrxnname(sim,rname,&order,&rxn,NULL,3);
1782 CHECKS(rxn,"unrecognized reaction name");
1783 CHECKS(order==2,"reaction_intersurface is only for bimolecular reactions");
1784 if(rxn->nprod==0) {
1785 line2=strnword(line2,2);
1786 CHECKS(line2,"no value entered");
1787 itct=sscanf(line2,"%s",nm);
1788 CHECKS(itct==1,"unable to read value");
1789 if(!strcmp(nm,"on")) rulelist[0]=0;
1790 else if(!strcmp(nm,"off")) rulelist[0]=-1;
1791 else CHECKS(0,"value not recognized"); }
1792 else {
1793 for(prd=0;prd<rxn->nprod;prd++) {
1794 line2=strnword(line2,2);
1795 CHECKS(line2,"not enough product codes entered");
1796 itct=sscanf(line2,"%s",nm);
1797 CHECKS(itct==1,"unable to read product code");
1798 if(prd==0 && !strcmp(nm,"off")) {
1799 rulelist[0]=-1;
1800 break; }
1801 if(!strcmp(nm,"+")) prd--;
1802 else if(nm[0]=='r') {
1803 if(nm[1]=='1') rulelist[prd]=1;
1804 else if(nm[1]=='2') rulelist[prd]=2;
1805 else CHECKS(0,"error reading a reactant code"); }
1806 else
1807 CHECKS(0,"error reading product rule"); }}
1808
1809 rxn=NULL;
1810 r=readrxnname(sim,rname,NULL,&rxn,NULL,1);
1811 if(r>=0) {
1812 er=RxnSetIntersurfaceRules(rxn,rulelist);
1813 CHECKS(er!=1,"out of memory allocating reaction intersurface rules"); }
1814 r=readrxnname(sim,rname,NULL,&rxn,&vlist,2);
1815 CHECKS(r!=-2,"out of memory reading reaction name");
1816 if(r>=0) {
1817 for(i=0;i<vlist->n;i++) {
1818 er=RxnSetIntersurfaceRules((rxnptr) vlist->xs[i],rulelist);
1819 CHECKS(er!=1,"out of memory allocating reaction intersurface rules"); }
1820 ListFreeV(vlist); }
1821 r=readrxnname(sim,rname,NULL,&rxn,NULL,3);
1822 CHECKS(r!=-2,"out of memory reading reaction name");
1823 if(r>=0) {
1824 er=RxnSetIntersurfaceRules(rxn,rulelist);
1825 CHECKS(er!=1,"out of memory allocating reaction intersurface rules"); }
1826 CHECKS(!strnword(line2,2),"unexpected text following reaction_intersurface"); }
1827
1828 else if(!strcmp(word,"reaction_representation")) { // reaction_representation
1829 itct=sscanf(line2,"%s",rname);
1830 CHECKS(itct==1,"reaction_representation format: rname rules_list");
1831 rxn=NULL;
1832 readrxnname(sim,rname,&order,&rxn,NULL,1);
1833 if(!rxn) readrxnname(sim,rname,&order,&rxn,NULL,2);
1834 if(!rxn) readrxnname(sim,rname,&order,&rxn,NULL,3);
1835 CHECKS(rxn,"unrecognized reaction name");
1836 for(rct=0;rct<order;rct++) {
1837 line2=strnword(line2,2);
1838 CHECKS(line2,"not enough reactant codes entered");
1839 itct=sscanf(line2,"%s",nm);
1840 CHECKS(itct==1,"unable to read reactant code");
1841 if(!strcmp(nm,"+")) rct--;
1842 else {
1843 replist[rct]=rxnstring2sr(nm);
1844 CHECKS(replist[rct]!=SRnone,"error reading a reactant code"); }}
1845 line2=strnword(line2,2);
1846 CHECKS(line2,"missing reaction arrow");
1847 itct=sscanf(line2,"%s",nm);
1848 CHECKS(itct==1,"unable to find reaction arrow");
1849 CHECKS(!strcmp(nm,"->"),"unable to find reaction arrow");
1850 for(prd=0;prd<rxn->nprod;prd++) {
1851 line2=strnword(line2,2);
1852 CHECKS(line2,"not enough product codes entered");
1853 itct=sscanf(line2,"%s",nm);
1854 CHECKS(itct==1,"unable to read product code");
1855 if(!strcmp(nm,"+")) prd--;
1856 else {
1857 replist[order+prd]=rxnstring2sr(nm);
1858 CHECKS(replist[order+prd]!=SRboth,"product representation cannot be 'both'");
1859 CHECKS(replist[order+prd]!=SRnone,"error reading a product code"); }}
1860
1861 rxn=NULL;
1862 r=readrxnname(sim,rname,&order,&rxn,NULL,1);
1863 if(r>=0) {
1864 er=RxnSetRepresentationRules(rxn,order,replist,replist+order);
1865 CHECKS(er!=1,"out of memory allocating reaction representation rules"); }
1866 r=readrxnname(sim,rname,&order,&rxn,&vlist,2);
1867 CHECKS(r!=-2,"out of memory reading reaction name");
1868 if(r>=0) {
1869 for(i=0;i<vlist->n;i++) {
1870 er=RxnSetRepresentationRules((rxnptr) vlist->xs[i],order,replist,replist+order);
1871 CHECKS(er!=1,"out of memory allocating reaction representation rules"); }
1872 ListFreeV(vlist); }
1873 r=readrxnname(sim,rname,&order,&rxn,NULL,3);
1874 CHECKS(r!=-2,"out of memory reading reaction name");
1875 if(r>=0) {
1876 er=RxnSetRepresentationRules(rxn,order,replist,replist+order);
1877 CHECKS(er!=1,"out of memory allocating reaction representation rules"); }
1878 CHECKS(!strnword(line2,2),"unexpected text following reaction_representation"); }
1879
1880 else if(!strcmp(word,"reaction_log")) { // reaction_log
1881 itct=sscanf(line2,"%s %s %s",fname,rname,nm);
1882 CHECKS(itct==3,"reaction_log format: filename rxnname serial_numbers");
1883 lilist=NULL;
1884 if(!strcmp(nm,"all")) {
1885 lilist=ListAppendItemLI(NULL,-1);
1886 CHECKS(lilist,"out of memory in reaction_log");
1887 CHECKS(!strnword(line2,4),"unexpected text following reaction_log"); }
1888 else {
1889 line2=strnword(line2,3);
1890 lilist=ListReadStringLI(line2);
1891 CHECKS(lilist,"unable to read serial numbers"); }
1892
1893 if(!strcmp(rname,"all")) { // all current reactions, no rules
1894 er=RxnSetLog(sim,fname,NULL,lilist,1);
1895 CHECKS(er!=2,"cannot log a reaction to multiple output files");
1896 CHECKS(!er,"out or memory in reaction_log"); }
1897 else {
1898 rxn=NULL;
1899 r=readrxnname(sim,rname,NULL,&rxn,NULL,1);
1900 if(r>=0) {
1901 er=RxnSetLog(sim,fname,rxn,lilist,1);
1902 CHECKS(er!=2,"cannot log a reaction to multiple output files");
1903 CHECKS(!er,"out or memory in reaction_log"); }
1904 r=readrxnname(sim,rname,NULL,&rxn,&vlist,2);
1905 CHECKS(r!=-2,"out of memory reading reaction name");
1906 if(r>=0) {
1907 for(i=0;i<vlist->n;i++) {
1908 er=RxnSetLog(sim,fname,(rxnptr) vlist->xs[i],lilist,1);
1909 CHECKS(er!=2,"cannot log a reaction to multiple output files");
1910 CHECKS(!er,"out or memory in reaction_log"); }
1911 ListFreeV(vlist); }
1912 r=readrxnname(sim,rname,NULL,&rxn,NULL,3);
1913 CHECKS(r!=-2,"out of memory reading reaction name");
1914 if(r>=0) {
1915 er=RxnSetLog(NULL,fname,rxn,lilist,1);
1916 CHECKS(er!=2,"cannot log a reaction to multiple output files");
1917 CHECKS(!er,"out or memory in reaction_log"); }}
1918 ListFreeLI(lilist);
1919 lilist=NULL; }
1920
1921 else if(!strcmp(word,"reaction_log_off")) { // reaction_log_off
1922 itct=sscanf(line2,"%s %s",rname,nm);
1923 CHECKS(itct==2,"reaction_log_off format: rxnname serial_numbers");
1924 lilist=NULL;
1925 if(!strcmp(nm,"all")) {
1926 lilist=ListAppendItemLI(NULL,-1);
1927 CHECKS(lilist,"out of memory in reaction_log");
1928 CHECKS(!strnword(line2,3),"unexpected text following reaction_log_off"); }
1929 else {
1930 line2=strnword(line2,2);
1931 lilist=ListReadStringLI(line2);
1932 CHECKS(lilist,"unable to read serial numbers"); }
1933
1934 if(!strcmp(rname,"all")) {
1935 er=RxnSetLog(sim,fname,NULL,lilist,0);
1936 CHECKS(!er,"out or memory in reaction_log"); }
1937 else {
1938 rxn=NULL;
1939 r=readrxnname(sim,rname,NULL,&rxn,NULL,1);
1940 if(r>=0) {
1941 er=RxnSetLog(sim,NULL,rxn,lilist,0);
1942 CHECKS(!er,"out or memory in reaction_log_off"); }
1943 r=readrxnname(sim,rname,NULL,&rxn,&vlist,2);
1944 CHECKS(r!=-2,"out of memory reading reaction name");
1945 if(r>=0) {
1946 for(i=0;i<vlist->n;i++) {
1947 er=RxnSetLog(sim,NULL,(rxnptr) vlist->xs[i],lilist,0);
1948 CHECKS(!er,"out or memory in reaction_log_off"); }
1949 ListFreeV(vlist); }
1950 r=readrxnname(sim,rname,NULL,&rxn,NULL,3);
1951 CHECKS(r!=-2,"out of memory reading reaction name");
1952 if(r>=0) {
1953 er=RxnSetLog(NULL,NULL,rxn,lilist,0);
1954 CHECKS(!er,"out or memory in reaction_log_off"); }}
1955 ListFreeLI(lilist);
1956 lilist=NULL; }
1957
1958 else if(!strcmp(word,"product_placement")) { // product_placement
1959 itct=sscanf(line2,"%s %s",rname,nm1);
1960 CHECKS(itct==2,"product_placement format: rname type parameters");
1961 rxn=NULL;
1962 readrxnname(sim,rname,&order,&rxn,NULL,1);
1963 if(!rxn) readrxnname(sim,rname,&order,&rxn,NULL,2);
1964 if(!rxn) readrxnname(sim,rname,&order,&rxn,NULL,3);
1965 CHECKS(rxn,"unrecognized reaction name");
1966 rpart=rxnstring2rp(nm1);
1967 CHECKS(!(rpart==RPbounce && (order!=2 || rxn->nprod<2)),"bounce can only be used with two reactants and at least two products");
1968 flt1=0;
1969 prd=0;
1970 for(d=0;d<sim->dim;d++) v1[prd]=0;
1971 line2=strnword(line2,3);
1972 if(rpart==RPirrev);
1973 else if(rpart==RPbounce && !line2)
1974 flt1=-2;
1975 else if(rpart==RPpgem || rpart==RPbounce || rpart==RPpgemmax || rpart==RPpgemmaxw || rpart==RPratio || rpart==RPunbindrad || rpart==RPpgem2 || rpart==RPpgemmax2 || rpart==RPratio2) {
1976 CHECKS(line2,"missing parameter in product_placement");
1977 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
1978 CHECKS(itct==1,"error reading parameter in product_placement");
1979 line2=strnword(line2,2); }
1980 else if(rpart==RPoffset || rpart==RPfixed) {
1981 CHECKS(line2,"missing parameters in product_placement");
1982 itct=sscanf(line2,"%s",nm1);
1983 CHECKS(itct==1,"format for product_placement: rname type parameters");
1984 if(strbegin("product_",nm1,0)) {
1985 itct=strmathsscanf(nm1+8,"%mi",varnames,varvalues,nvar,&prd);
1986 CHECKS(itct==1,"error reading product number");
1987 prd--;
1988 CHECKS(prd>=0 && prd<rxn->nprod,"illegal product number"); }
1989 else {
1990 i=molstring2index1(sim,nm1,&ms,&index);
1991 CHECKS(i!=0,"wildcards are not permitted in product_placement");
1992 CHECKS(i!=-1,"cannot read molecule species name");
1993 CHECKS(i!=-2,"mismatched or improper parentheses around product state");
1994 CHECKS(i!=-3,"cannot read molecule state value");
1995 CHECKS(i!=-4,"molecule name not recognized");
1996 CHECKS(i!=-5,"'all' is not allowed in product_placement");
1997 CHECKS(i!=-6,"cannot read wildcard logic");
1998 CHECKS(i!=-7,"error allocating memory");
1999 CHECKS(ms<MSMAX1 || ms==MSall,"invalid state");
2000 for(prd=0;prd<rxn->nprod && rxn->prdident[prd]!=i && rxn->prdstate[prd]!=ms;prd++);
2001 CHECKS(prd<rxn->nprod,"molecule in product_placement is not a product of this reaction"); }
2002 CHECKS(line2=strnword(line2,2),"position vector missing for product_placement");
2003 for(d=0;d<dim;d++) {
2004 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&v1[d]);
2005 CHECKS(itct==1,"insufficient data for position vector for product_placement");
2006 line2=strnword(line2,2); }}
2007 else CHECKS(0,"unrecognized or not permitted product placement parameter");
2008
2009 rxn=NULL;
2010 r=readrxnname(sim,rname,NULL,&rxn,NULL,1);
2011 if(r>=0) {
2012 i1=RxnSetRevparam(sim,rxn,rpart,flt1,prd,v1,dim);
2013 CHECKS(i1!=2,"reversible parameter value is out of bounds"); }
2014 r=readrxnname(sim,rname,NULL,&rxn,&vlist,2);
2015 CHECKS(r!=-2,"out of memory reading reaction name");
2016 if(r>=0) {
2017 for(i=0;i<vlist->n;i++) {
2018 i1=RxnSetRevparam(sim,(rxnptr) vlist->xs[i],rpart,flt1,prd,v1,dim);
2019 CHECKS(i1!=2,"reversible parameter value is out of bounds"); }
2020 ListFreeV(vlist); }
2021 r=readrxnname(sim,rname,NULL,&rxn,NULL,3);
2022 CHECKS(r!=-2,"out of memory reading reaction name");
2023 if(r>=0) {
2024 i1=RxnSetRevparam(NULL,rxn,rpart,flt1,prd,v1,dim);
2025 CHECKS(i1!=2,"reversible parameter value is out of bounds"); }
2026 CHECKS(!line2,"unexpected text following product_placement"); }
2027
2028 else if(!strcmp(word,"expand_rules")) { // expand_rules
2029 itct=sscanf(line2,"%s",nm);
2030 if(!strcmp(nm,"all")) i=-1;
2031 else if(!strcmp(nm,"otf") || !strcmp(nm,"on-the-fly")) i=-2;
2032 else {
2033 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&i);
2034 CHECKS(itct==1,"expand_rules format: iterations");
2035 CHECKS(i>=0,"iterations needs to be >=0"); }
2036 er=RuleExpandRules(sim,i);
2037 CHECKS(er!=-1 && er!=-11,"out of memory while expanding rules");
2038 CHECKS(er!=-5,"generated species name exceeded maximum allowed length of %i characters",STRCHAR);
2039 CHECKS(er!=-41,"no rules to expand");
2040 CHECKS(!er,"BUG: error number %i arose in RuleExpandRules function",er);
2041 CHECKS(!strnword(line2,2),"unexpected text following expand_rules"); }
2042
2043 // optimizing runtime
2044
2045 else if(!strcmp(word,"rand_seed") || !strcmp(word,"random_seed")) { // rand_seed, random_seed
2046 itct=sscanf(line2,"%s",nm);
2047 if(!strcmp(nm,"time"))
2048 li1=(long int) time(NULL);
2049 else {
2050 itct=sscanf(line2,"%li",&li1);
2051 CHECKS(itct==1,"random_seed needs to be 'time' or an integer"); }
2052 Simsetrandseed(sim,li1);
2053 CHECKS(!strnword(line2,2),"unexpected text following random_seed"); }
2054
2055 else if(!strcmp(word,"accuracy")) { // accuracy
2056 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
2057 CHECKS(itct==1,"accuracy needs to be a number");
2058 sim->accur=flt1;
2059 CHECKS(!strnword(line2,2),"unexpected text following accuracy"); }
2060
2061 else if(!strcmp(word,"molperbox")) { // molperbox
2062 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
2063 CHECKS(itct==1,"molperbox needs to be a number");
2064 er=boxsetsize(sim,"molperbox",flt1);
2065 CHECKS(er!=1,"out of memory");
2066 CHECKS(er!=2,"molperbox needs to be >0");
2067 CHECKS(er!=3,"need to enter dim before molperbox");
2068 CHECKS(!strnword(line2,2),"unexpected text following molperbox"); }
2069
2070 else if(!strcmp(word,"boxsize")) { // boxsize
2071 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
2072 CHECKS(itct==1,"boxsize needs to be a number");
2073 er=boxsetsize(sim,"boxsize",flt1);
2074 CHECKS(er!=1,"out of memory");
2075 CHECKS(er!=2,"boxsize needs to be >0");
2076 CHECKS(er!=3,"need to enter dim before boxsize");
2077 CHECKS(!strnword(line2,2),"unexpected text following boxsize"); }
2078
2079 else if(!strcmp(word,"gauss_table_size")) { // gauss_table_size
2080 itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&i1);
2081 CHECKS(itct==1,"gauss_table_size needs to be an integer");
2082 er=molssetgausstable(sim,i1);
2083 CHECKS(er!=1,"out of memory");
2084 CHECKS(er!=2,"need to enter max_species before gauss_table_size");
2085 CHECKS(er!=3,"gauss_table_size needs to be an integer power of two");
2086 CHECKS(!strnword(line2,2),"unexpected text following gauss_table_size"); }
2087
2088 else if(!strcmp(word,"epsilon")) { // epsilon
2089 CHECKS(dim>0,"need to enter dim before epsilon");
2090 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
2091 CHECKS(itct==1,"epsilon format: value");
2092 er=surfsetepsilon(sim,flt1);
2093 CHECKS(er!=2,"out of memory");
2094 CHECKS(er!=3,"epsilon value needs to be at least 0");
2095 CHECKS(!strnword(line2,2),"unexpected text following epsilon"); }
2096
2097 else if(!strcmp(word,"margin")) { // margin
2098 CHECKS(dim>0,"need to enter dim before margin");
2099 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
2100 CHECKS(itct==1,"margin format: value");
2101 er=surfsetmargin(sim,flt1);
2102 CHECKS(er!=2,"out of memory");
2103 CHECKS(er!=3,"margin value needs to be at least 0");
2104 CHECKS(!strnword(line2,2),"unexpected text following margin"); }
2105
2106 else if(!strcmp(word,"neighbor_dist") || !strcmp(word,"neighbour_dist")) { // neighbor_dist
2107 CHECKS(dim>0,"need to enter dim before neighbor_dist");
2108 itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&flt1);
2109 CHECKS(itct==1,"neighbor_dist format: value");
2110 er=surfsetneighdist(sim,flt1);
2111 CHECKS(er!=2,"out of memory");
2112 CHECKS(er!=3,"neighdist value needs to be at least 0");
2113 CHECKS(!strnword(line2,2),"unexpected text following neighbor_dist"); }
2114
2115 else { // unknown word
2116 CHECKS(0,"syntax error: statement not recognized"); }
2117
2118 return 0;
2119
2120 failure:
2121 simParseError(sim,pfp);
2122 return 1; }
2123
2124
2125 #ifdef VCELL
2126 int loadJMS(simptr sim,ParseFilePtr *pfpptr,char *line2,char *erstr);
2127 #endif
2128
2129
2130 /* loadsim */
loadsim(simptr sim,const char * fileroot,const char * filename,const char * flags)2131 int loadsim(simptr sim,const char *fileroot,const char *filename,const char *flags) {
2132 int done,pfpcode,er;
2133 char word[STRCHAR],*line2,errstring[STRCHARLONG];
2134 ParseFilePtr pfp;
2135
2136 strncpy(sim->filepath,fileroot,STRCHAR);
2137 strncpy(sim->filename,filename,STRCHAR);
2138 strncpy(sim->flags,flags,STRCHAR);
2139 strcpy(SimFlags,flags);
2140 done=0;
2141 ErrorLineAndString[0]='\0';
2142
2143 pfp=Parse_Start(fileroot,filename,errstring);
2144 CHECKS(pfp,"%s",errstring);
2145 er=Parse_CmdLineArg(NULL,NULL,pfp);
2146 CHECKMEM(!er);
2147 #ifdef OPTION_VCELL
2148 sim->volumeSamplesPtr = NULL;//initialize the volumesample to null
2149 #endif
2150
2151 while(!done) {
2152 if(pfp->lctr==0 && !strchr(flags,'q'))
2153 simLog(sim,2," Reading file: '%s'\n",pfp->fname);
2154 pfpcode=Parse_ReadLine(&pfp,word,&line2,errstring);
2155 CHECKS(pfpcode!=3,"%s",errstring);
2156 er=0;
2157
2158 if(pfpcode==0); // already taken care of
2159
2160 else if(pfpcode==2) { // end reading
2161 done=1; }
2162
2163 #ifdef VCELL
2164 else if(!strcmp(word,"start_jms")) { // jms settings
2165 er=loadJMS(sim,&pfp,line2, errstring);
2166 }
2167
2168 else if(!strcmp(word,"start_highResVolumeSamples")) { //highResVolumeSamplesFile
2169 er=loadHighResVolumeSamples(sim,&pfp,line2); }
2170 #endif
2171
2172 else if(!strcmp(word,"start_reaction")) { // start_reaction
2173 CHECKS(1,"reactions cannot be entered in blocks; this feature was deprecated"); }
2174
2175 else if(!strcmp(word,"start_surface")) { // start_surface
2176 CHECKS(sim->dim>0,"need to enter dim before start_surface");
2177 er=loadsurface(sim,&pfp,line2); }
2178
2179 else if(!strcmp(word,"start_compartment")) { // start_compartment
2180 CHECKS(sim->dim>0,"need to enter dim before start_compartment");
2181 er=loadcompart(sim,&pfp,line2); }
2182
2183 else if(!strcmp(word,"start_port")) { // start_port
2184 er=loadport(sim,&pfp,line2); }
2185
2186 else if(!strcmp(word,"start_lattice")) { // start_lattice
2187 er=loadlattice(sim,&pfp,line2); }
2188
2189 else if(!strcmp(word,"start_bng")) { // start_bng
2190 er=loadbng(sim,&pfp,line2);
2191 if(er) return 1;
2192 er=bngupdate(sim); }
2193
2194 else if(!strcmp(word,"start_filament_type")) { // start_filament_type
2195 er=filloadtype(sim,&pfp,line2); }
2196
2197 else if(!strcmp(word,"start_filament")) { // start_filament
2198 er=filloadfil(sim,&pfp,line2,NULL); }
2199
2200 else if(!strcmp(word,"start_rules")) { // start_rules
2201 CHECKS(0,"Moleculizer support has been discontinued in Smoldyn"); }
2202
2203 else if(!line2) { // just word
2204 CHECKS(0,"unknown word or missing parameter"); }
2205
2206 else {
2207 er=simreadstring(sim,pfp,word,line2); }
2208
2209 if(er) return 1; }
2210
2211 return 0;
2212
2213 failure: // failure
2214 simParseError(sim,pfp);
2215 return 1; }
2216
2217
2218 /* simupdate */
simupdate(simptr sim)2219 int simupdate(simptr sim) {
2220 int er;
2221 static int recurs=0;
2222
2223 if(sim->condition==SCok) {
2224 return 0; }
2225 if(recurs>10) {
2226 recurs=0;
2227 return 2; }
2228 recurs++;
2229
2230 if(sim->condition==SCinit && sim->mols)
2231 simLog(sim,2," setting up molecules\n");
2232 er=molsupdate(sim);
2233 CHECK(er!=1);
2234
2235 if(sim->condition==SCinit)
2236 simLog(sim,2," setting up virtual boxes\n");
2237 er=boxesupdate(sim);
2238 CHECK(er!=1);
2239 CHECKS(er!=3,"simulation dimensions or boundaries are undefined");
2240
2241 er=molsort(sim,0);
2242 CHECK(er!=1);
2243
2244 if(sim->condition==SCinit && sim->cmptss)
2245 simLog(sim,2," setting up compartments\n");
2246 er=compartsupdate(sim);
2247 CHECK(er!=1);
2248
2249 if(sim->condition==SCinit && (sim->rxnss[0] || sim->rxnss[1] || sim->rxnss[2]))
2250 simLog(sim,2," setting up reactions\n");
2251 er=rxnsupdate(sim);
2252 CHECK(er!=1);
2253 CHECKS(er!=3,"failed to set up reactions");
2254
2255 if(sim->condition==SCinit && sim->srfss)
2256 simLog(sim,2," setting up surfaces\n");
2257 er=surfupdate(sim);
2258 CHECK(er!=1);
2259
2260 if(sim->condition==SCinit && sim->portss)
2261 simLog(sim,2," setting up ports\n");
2262 er=portsupdate(sim);
2263 CHECK(er!=1);
2264
2265 if(sim->condition==SCinit && sim->latticess)
2266 simLog(sim,2," setting up lattices\n");
2267 er=latticesupdate(sim);
2268 CHECK(er!=1);
2269
2270 if(sim->condition==SCinit && sim->filss)
2271 simLog(sim,2," setting up filaments\n");
2272 er=filsupdate(sim);
2273 CHECK(er!=1);
2274
2275 if(sim->condition==SCinit && sim->graphss)
2276 simLog(sim,2," setting up graphics\n");
2277 er=graphicsupdate(sim);
2278 CHECK(er!=1);
2279
2280 if(sim->mols && sim->mols->condition!=SCok) {
2281 er=simupdate(sim);
2282 CHECK(!er); }
2283 if(sim->boxs && sim->boxs->condition!=SCok) {
2284 er=simupdate(sim);
2285 CHECK(!er); }
2286 if(sim->cmptss && sim->cmptss->condition!=SCok) {
2287 er=simupdate(sim);
2288 CHECK(!er); }
2289 if(sim->rxnss[0] && sim->rxnss[0]->condition!=SCok) {
2290 er=simupdate(sim);
2291 CHECK(!er); }
2292 if(sim->rxnss[1] && sim->rxnss[1]->condition!=SCok) {
2293 er=simupdate(sim);
2294 CHECK(!er); }
2295 if(sim->rxnss[2] && sim->rxnss[2]->condition!=SCok) {
2296 er=simupdate(sim);
2297 CHECK(!er); }
2298 if(sim->srfss && sim->srfss->condition!=SCok) {
2299 er=simupdate(sim);
2300 CHECK(!er); }
2301 if(sim->portss && sim->portss->condition!=SCok) {
2302 er=simupdate(sim);
2303 CHECK(!er); }
2304 if(sim->filss && sim->filss->condition!=SCok) {
2305 er=simupdate(sim);
2306 CHECK(!er); }
2307 if(sim->graphss && sim->graphss->condition!=SCok) {
2308 er=simupdate(sim);
2309 CHECK(!er); }
2310
2311 er=reassignmolecs(sim,0,0);
2312 CHECK(!er);
2313
2314 if(sim->cmds && sim->cmds->condition!=3) {
2315 er=scmdupdatecommands(sim->cmds,sim->tmin,sim->tmax,sim->dt);
2316 CHECK(!er); }
2317
2318 simsetcondition(sim,SCok,1);
2319 recurs=0;
2320
2321 return 0;
2322
2323 failure:
2324 if(ErrorType!=1) simLog(sim,10,"%s",ErrorString);
2325 return 1; }
2326
2327
2328 /* simInitAndLoad */
2329 #ifdef OPTION_VCELL
simInitAndLoad(const char * fileroot,const char * filename,simptr * smptr,const char * flags,ValueProviderFactory * valueProviderFactory,AbstractMesh * mesh)2330 int simInitAndLoad(const char *fileroot,const char *filename,simptr *smptr,const char *flags, ValueProviderFactory* valueProviderFactory, AbstractMesh* mesh) {
2331
2332 simptr sim;
2333 int er,qflag,sflag;
2334
2335 sim=*smptr;
2336 if(!sim) {
2337 qflag=strchr(flags,'q')?1:0;
2338 sflag=strchr(flags,'s')?1:0;
2339 if(!qflag && !sflag) {
2340 simLog(NULL,2,"--------------------------------------------------------------\n");
2341 simLog(NULL,2,"Running Smoldyn %s\n",VERSION);
2342 simLog(NULL,2,"\nCONFIGURATION FILE\n");
2343 simLog(NULL,2," Path: '%s'\n",fileroot);
2344 simLog(NULL,2," Name: '%s'\n",filename); }
2345 sim=simalloc(fileroot);
2346 CHECKMEM(sim);
2347 sim->valueProviderFactory = valueProviderFactory; //create value provider factory
2348 sim->mesh = mesh;
2349 er=loadsim(sim,fileroot,filename,flags); // load sim
2350 CHECK(!er);
2351 if(sim && sim->valueProviderFactory) {
2352 sim->valueProviderFactory->setSimptr(sim); }
2353
2354 simLog(sim,2," Loaded file successfully\n"); }
2355
2356 *smptr=sim;
2357 return 0;
2358
2359 failure:
2360 if(ErrorType!=1) simLog(sim,10,ErrorString);
2361 if(!*smptr) simfree(sim);
2362 return 1; }
2363
2364 #else
2365
2366
2367 /* simInitAndLoad */
simInitAndLoad(const char * fileroot,const char * filename,simptr * smptr,const char * flags)2368 int simInitAndLoad(const char *fileroot,const char *filename,simptr *smptr,const char *flags) {
2369 simptr sim;
2370 int er,qflag,sflag;
2371
2372 sim=*smptr;
2373 if(!sim) {
2374 qflag=strchr(flags,'q')?1:0;
2375 sflag=strchr(flags,'s')?1:0;
2376 if(!qflag && !sflag) {
2377 simLog(NULL,2,"--------------------------------------------------------------\n");
2378 simLog(NULL,2,"Running Smoldyn %s\n",VERSION);
2379 simLog(NULL,2,"\nCONFIGURATION FILE\n");
2380 simLog(NULL,2," Path: '%s'\n",fileroot);
2381 simLog(NULL,2," Name: '%s'\n",filename); }
2382 sim=simalloc(fileroot);
2383 CHECKMEM(sim);
2384 CHECKMEM(strloadmathfunctions()==0);
2385 CHECKMEM(loadsmolfunctions(sim)==0);
2386 er=loadsim(sim,fileroot,filename,flags); // load sim
2387 CHECK(!er);
2388 simLog(sim,2," Loaded file successfully\n"); }
2389
2390 *smptr=sim;
2391 return 0;
2392
2393 failure:
2394 if(ErrorType!=1) simLog(sim,10,ErrorString);
2395 if(!*smptr) simfree(sim);
2396 return 1; }
2397
2398 #endif
2399
2400
2401 /* simUpdateAndDisplay */
simUpdateAndDisplay(simptr sim)2402 int simUpdateAndDisplay(simptr sim) {
2403 int er;
2404
2405 er=simupdate(sim); // update sim
2406 CHECK(!er);
2407
2408 simLog(sim,2,"\n"); // diagnostics output
2409 simsystemoutput(sim);
2410 er=checksimparams(sim);
2411 CHECK(!er);
2412 return 0;
2413
2414 failure:
2415 return 1; }
2416
2417
2418 /******************************************************************************/
2419 /************************** core simulation functions *************************/
2420 /******************************************************************************/
2421
2422
2423 /* simdocommands */
simdocommands(simptr sim)2424 int simdocommands(simptr sim) {
2425 int er;
2426 enum CMDcode ccode;
2427
2428 ccode=scmdexecute(sim->cmds,sim->time,sim->dt,-1,0);
2429 er=simupdate(sim);
2430 if(er) return 8;
2431 er=molsort(sim,0); // sort live and dead
2432 if(er) return 6;
2433 if(ccode==CMDstop || ccode==CMDabort) return 7;
2434 return 0; }
2435
2436
2437 /* debugcode */
debugcode(simptr sim,const char * prefix)2438 void debugcode(simptr sim,const char *prefix) {
2439 int m;
2440 moleculeptr mptr;
2441 char string[STRCHAR];
2442
2443 if(sim->time<189.243 || sim->time>189.247) return;
2444 for(m=0;m<sim->mols->nl[0];m++) {
2445 mptr=sim->mols->live[0][m];
2446 if(mptr->serno==1377166 || mptr->serno==1374858) {
2447 printf("%s: time=%g serno=%s",prefix,sim->time,molserno2string(mptr->serno,string));
2448 printf(" posx=(%g,%g,%g)",mptr->posx[0],mptr->posx[1],mptr->posx[2]);
2449 printf(" pos=(%g,%g,%g)",mptr->pos[0],mptr->pos[1],mptr->pos[2]);
2450 printf(" pnl=%s",mptr->pnl?mptr->pnl->pname:"NULL");
2451 printf(" posx side=%s",surfface2string(panelside(mptr->posx,sim->srfss->srflist[4]->panels[PSdisk][0],3,NULL,1,0),string));
2452 printf(" pos side=%s",surfface2string(panelside(mptr->pos,sim->srfss->srflist[4]->panels[PSdisk][0],3,NULL,1,0),string));
2453 printf("\n"); }}
2454 return; }
2455
2456
2457
2458 /* simulatetimestep */
simulatetimestep(simptr sim)2459 int simulatetimestep(simptr sim) {
2460 int er,ll;
2461
2462 er=RuleExpandRules(sim,-3); // expand any reaction rules if needed
2463 if(er && er!=-41) return 13;
2464
2465 er=simupdate(sim); // update any data structure changes
2466 if(er) return 8;
2467
2468 er=(*sim->diffusefn)(sim); // diffuse
2469 if(er) return 9;
2470
2471 if(sim->srfss) { // deal with surface or wall collisions
2472 for(ll=0;ll<sim->srfss->nmollist;ll++) {
2473 if(sim->srfss->srfmollist[ll] & SMLdiffuse) {
2474 (*sim->surfacecollisionsfn)(sim,ll,0); }}
2475 for(ll=0;ll<sim->srfss->nmollist;ll++) // surface-bound molecule actions
2476 if(sim->srfss->srfmollist[ll] & SMLsrfbound)
2477 (*sim->surfaceboundfn)(sim,ll); }
2478 else {
2479 if(sim->mols)
2480 for(ll=0;ll<sim->mols->nlist;ll++)
2481 if(sim->mols->diffuselist[ll])
2482 (*sim->checkwallsfn)(sim,ll,0,NULL); }
2483
2484 er=(*sim->assignmols2boxesfn)(sim,1,0); // assign to boxes (diffusing molecs., not reborn)
2485 if(er) return 2;
2486
2487 er=molsort(sim,0); // sort live and dead
2488 if(er) return 6;
2489
2490 er=(*sim->zeroreactfn)(sim);
2491 if(er) return 3;
2492
2493 er=(*sim->unimolreactfn)(sim);
2494 if(er) return 4;
2495
2496 er=(*sim->bimolreactfn)(sim,0);
2497 if(er) return 5;
2498 er=(*sim->bimolreactfn)(sim,1);
2499 if(er) return 5;
2500
2501 er=molsort(sim,0); // sort live and dead
2502 if(er) return 6;
2503
2504 if(sim->latticess) {
2505 er=latticeruntimestep(sim);
2506 if(er) return 12;
2507 er=molsort(sim,1);
2508 if(er) return 6; }
2509
2510 if(sim->srfss) {
2511 for(ll=0;ll<sim->srfss->nmollist;ll++)
2512 (*sim->surfacecollisionsfn)(sim,ll,1); } // surfaces again, reborn molecs. only
2513 else {
2514 if(sim->mols)
2515 for(ll=0;ll<sim->mols->nlist;ll++)
2516 (*sim->checkwallsfn)(sim,ll,1,NULL); }
2517
2518 er=(*sim->assignmols2boxesfn)(sim,0,1); // assign again (all, reborn)
2519 if(er) return 2;
2520
2521 er=filDynamics(sim);
2522 if(er) return 11;
2523
2524 #ifdef ENABLE_PYTHON_CALLBACK
2525 for(unsigned int i = 0; i < sim->ncallbacks; i++) {
2526 auto callback = sim->callbacks[i];
2527 if((0 == (sim->simstep % callback->getStep())) && callback->isValid()) {
2528 callback->evalAndUpdate(sim->time);
2529 }
2530 }
2531 sim->simstep += 1;
2532 #endif
2533 sim->time+=sim->dt; // --- end of time step ---
2534 simsetvariable(sim,"time",sim->time);
2535 er=simdocommands(sim);
2536 if(er) return er;
2537
2538
2539 if(sim->time>=sim->tmax) return 1;
2540 if(sim->time>=sim->tbreak) return 10;
2541
2542 return 0; }
2543
2544
2545 /* endsimulate */
endsimulate(simptr sim,int er)2546 void endsimulate(simptr sim,int er) {
2547 int sflag,tflag,*eventcount;
2548
2549 gl2State(2);
2550 tflag=strchr(sim->flags,'t')?1:0;
2551 sflag=strchr(sim->flags,'s')?1:0;
2552
2553 simLog(sim,2,"\n");
2554 if(er==1) simLog(sim,2,"Simulation complete\n");
2555 else if(er==2) simLog(sim,5,"Simulation terminated during molecule assignment\n Out of memory\n");
2556 else if(er==3) simLog(sim,5,"Simulation terminated during order 0 reaction\n Not enough molecules allocated\n Maximum allowed molecule number is %i",sim->mols->maxdlimit);
2557 else if(er==4) simLog(sim,5,"Simulation terminated during order 1 reaction\n Not enough molecules allocated\n Maximum allowed molecule number is %i",sim->mols->maxdlimit);
2558 else if(er==5) simLog(sim,5,"Simulation terminated during order 2 reaction\n Not enough molecules allocated\n Maximum allowed molecule number is %i",sim->mols->maxdlimit);
2559 else if(er==6) simLog(sim,5,"Simulation terminated during molecule sorting\n Out of memory\n");
2560 else if(er==7) simLog(sim,5,"Simulation stopped by a runtime command\n");
2561 else if(er==8) simLog(sim,5,"Simulation terminated during simulation state updating\n Out of memory\n");
2562 else if(er==9) simLog(sim,5,"Simulation terminated during diffusion\n Out of memory\n");
2563 else if(er==11) simLog(sim,5,"Simulation terminated during filament dynamics\n");
2564 else if(er==12) simLog(sim,5,"Simulation terminated during lattice simulation\n");
2565 else if(er==13) simLog(sim,5,"Simulation terminated during reaction network expansion\n");
2566 else simLog(sim,2,"Simulation stopped by user\n");
2567 simLog(sim,2,"Current simulation time: %f\n",sim->time);
2568
2569 eventcount=sim->eventcount;
2570 if(eventcount[ETwall]) simLog(sim,2,"%i wall interactions\n",eventcount[ETwall]);
2571 if(eventcount[ETsurf]) simLog(sim,2,"%i surface interactions\n",eventcount[ETsurf]);
2572 if(eventcount[ETdesorb]) simLog(sim,2,"%i desorptions\n",eventcount[ETdesorb]);
2573 if(eventcount[ETrxn0]) simLog(sim,2,"%i zeroth order reactions\n",eventcount[ETrxn0]);
2574 if(eventcount[ETrxn1]) simLog(sim,2,"%i unimolecular reactions\n",eventcount[ETrxn1]);
2575 if(eventcount[ETrxn2intra]) simLog(sim,2,"%i intrabox bimolecular reactions\n",eventcount[ETrxn2intra]);
2576 if(eventcount[ETrxn2inter]) simLog(sim,2,"%i interbox bimolecular reactions\n",eventcount[ETrxn2inter]);
2577 if(eventcount[ETrxn2wrap]) simLog(sim,2,"%i wrap-around bimolecular reactions\n",eventcount[ETrxn2wrap]);
2578 if(eventcount[ETrxn2hybrid]) simLog(sim,2,"%i bybrid bimolecular reactions\n",eventcount[ETrxn2hybrid]);
2579 if(eventcount[ETimport]) simLog(sim,2,"%i imported molecules\n",eventcount[ETimport]);
2580 if(eventcount[ETexport]) simLog(sim,2,"%i exported molecules\n",eventcount[ETexport]);
2581
2582 simLog(sim,2,"total execution time: %g seconds\n",sim->elapsedtime);
2583
2584 // TODO: Useful when running benchmarks
2585 // If SMOLDYN_NO_PROMPT is set by user then smoldyn quit at the
2586 // end of simultion without prompting user to press shift+Q.
2587 // It is useful when running tests in batch mode locally. And it is essential
2588 // for testing examples on remote servers such as Travis CI and github actions.
2589 const char* dontPrompt = getenv("SMOLDYN_NO_PROMPT");
2590 if(dontPrompt != NULL && strlen(dontPrompt) > 0)
2591 sim->quitatend = 1;
2592
2593 if(sim->graphss && sim->graphss->graphics>0 && !tflag && !sim->quitatend && !sflag)
2594 fprintf(stderr,"\nTo quit: Activate graphics window, then press shift-Q.\a\n");
2595 return; }
2596
2597
2598 /* smolsimulate */
smolsimulate(simptr sim)2599 int smolsimulate(simptr sim) {
2600 int er;
2601
2602 er=0;
2603 simLog(sim,2,"Simulating\n");
2604 sim->clockstt=time(NULL);
2605 er=simdocommands(sim);
2606 if(!er)
2607 while((er=simulatetimestep(sim))==0);
2608 if(er!=10) {
2609 scmdpop(sim->cmds,sim->tmax);
2610 scmdexecute(sim->cmds,sim->time,sim->dt,-1,1);
2611 scmdsetcondition(sim->cmds,0,0); }
2612 sim->elapsedtime+=difftime(time(NULL),sim->clockstt);
2613 return er; }
2614