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 <stdio.h>
9 #include <math.h>
10 #include <string.h>
11 #include "Geometry.h"
12 #include "List.h"
13 #include "math2.h"
14 #include "opengl2.h"
15 #include "random2.h"
16 #include "Rn.h"
17 #include "RnSort.h"
18 #include "string2.h"
19 #include "SimCommand.h"
20 #include "Zn.h"
21
22 #ifdef OPTION_NSV
23 #include "nsvc.h"
24 #endif
25
26 #include "smoldyn.h"
27 #include "smoldynfuncs.h"
28 #include "smoldynconfigure.h"
29
30
31 /* Global variables */
32 int Nvar;
33 char **Varnames;
34 double *Varvalues;
35
36
37 /**********************************************************/
38 /******************** command declarations ****************/
39 /**********************************************************/
40
41 #ifdef VCELL
42 #include <vcellcmd.h>
43 #endif
44
45 // simulation control
46 enum CMDcode cmdstop(simptr sim,cmdptr cmd,char *line2);
47 enum CMDcode cmdpause(simptr sim,cmdptr cmd,char *line2);
48 enum CMDcode cmdbeep(simptr sim,cmdptr cmd,char *line2);
49 enum CMDcode cmdkeypress(simptr sim,cmdptr cmd,char *line2);
50 enum CMDcode cmdsetflag(simptr sim,cmdptr cmd,char *line2);
51 enum CMDcode cmdsetrandseed(simptr sim,cmdptr cmd,char *line2);
52 enum CMDcode cmdsetgraphics(simptr sim,cmdptr cmd,char *line2);
53 enum CMDcode cmdsetgraphic_iter(simptr sim,cmdptr cmd,char *line2);
54 enum CMDcode cmdupdategraphics(simptr sim,cmdptr cmd,char *line2);
55
56 // file manipulation
57 enum CMDcode cmdoverwrite(simptr sim,cmdptr cmd,char *line2);
58 enum CMDcode cmdincrementfile(simptr sim,cmdptr cmd,char *line2);
59
60 // conditional
61 enum CMDcode cmdifflag(simptr sim,cmdptr cmd,char *line2);
62 enum CMDcode cmdifprob(simptr sim,cmdptr cmd,char *line2);
63 enum CMDcode cmdifno(simptr sim,cmdptr cmd,char *line2);
64 enum CMDcode cmdifless(simptr sim,cmdptr cmd,char *line2);
65 enum CMDcode cmdifmore(simptr sim,cmdptr cmd,char *line2);
66 enum CMDcode cmdifincmpt(simptr sim,cmdptr cmd,char *line2);
67 enum CMDcode cmdifchange(simptr sim,cmdptr cmd,char *line2);
68 enum CMDcode cmdif(simptr sim,cmdptr cmd,char *line2);
69
70 // system observation
71 enum CMDcode cmdecho(simptr sim,cmdptr cmd,char *line2);
72 enum CMDcode cmdevaluate(simptr sim,cmdptr cmd,char *line2);
73 enum CMDcode cmdwarnescapee(simptr sim,cmdptr cmd,char *line2);
74 enum CMDcode cmdwarnescapeecmpt(simptr sim,cmdptr cmd,char *line2);
75 enum CMDcode cmdmolcountheader(simptr sim,cmdptr cmd,char *line2);
76 enum CMDcode cmdmolcount(simptr sim,cmdptr cmd,char *line2);
77 enum CMDcode cmdmolcountinbox(simptr sim,cmdptr cmd,char *line2);
78 enum CMDcode cmdmolcountincmpt(simptr sim,cmdptr cmd,char *line2);
79 enum CMDcode cmdmolcountincmpts(simptr sim,cmdptr cmd,char *line2);
80 enum CMDcode cmdmolcountincmpt2(simptr sim,cmdptr cmd,char *line2);
81 enum CMDcode cmdmolcountonsurf(simptr sim,cmdptr cmd,char *line2);
82 enum CMDcode cmdmolcountspace(simptr sim,cmdptr cmd,char *line2);
83 enum CMDcode cmdmolcountspace2d(simptr sim,cmdptr cmd,char *line2);
84 enum CMDcode cmdmolcountspaceradial(simptr sim,cmdptr cmd,char *line2);
85 enum CMDcode cmdmolcountspacepolarangle(simptr sim,cmdptr cmd,char *line2);
86 enum CMDcode cmdradialdistribution(simptr sim,cmdptr cmd,char *line2);
87 enum CMDcode cmdradialdistribution2(simptr sim,cmdptr cmd,char *line2);
88 enum CMDcode cmdmolcountspecies(simptr sim,cmdptr cmd,char *line2);
89 enum CMDcode cmdmolcountspecieslist(simptr sim,cmdptr cmd,char *line2);
90 enum CMDcode cmdmollistsize(simptr sim,cmdptr cmd,char *line2);
91 enum CMDcode cmdlistmols(simptr sim,cmdptr cmd,char *line2);
92 enum CMDcode cmdlistmols2(simptr sim,cmdptr cmd,char *line2);
93 enum CMDcode cmdlistmols3(simptr sim,cmdptr cmd,char *line2);
94 enum CMDcode cmdlistmols4(simptr sim,cmdptr cmd,char *line2);
95 enum CMDcode cmdlistmolscmpt(simptr sim,cmdptr cmd,char *line2);
96 enum CMDcode cmdmolpos(simptr sim,cmdptr cmd,char *line2);
97 enum CMDcode cmdtrackmol(simptr sim,cmdptr cmd,char *line2);
98 enum CMDcode cmdmolmoments(simptr sim,cmdptr cmd,char *line2);
99 enum CMDcode cmdsavesim(simptr sim,cmdptr cmd,char *line2);
100 enum CMDcode cmdmeansqrdisp(simptr sim,cmdptr cmd,char *line2);
101 enum CMDcode cmdmeansqrdisp2(simptr sim,cmdptr cmd,char *line2);
102 enum CMDcode cmdmeansqrdisp3(simptr sim,cmdptr cmd,char *line2);
103 enum CMDcode cmdresidencetime(simptr sim,cmdptr cmd,char *line2);
104 enum CMDcode cmddiagnostics(simptr sim,cmdptr cmd,char *line2);
105 enum CMDcode cmdexecutiontime(simptr sim,cmdptr cmd,char *line2);
106 enum CMDcode cmdwriteVTK(simptr sim,cmdptr cmd,char *line2);
107 enum CMDcode cmdprintdata(simptr sim,cmdptr cmd,char *line2);
108
109 enum CMDcode cmdprintLattice(simptr sim,cmdptr cmd,char *line2);
110
111 // system manipulation
112 enum CMDcode cmdset(simptr sim,cmdptr cmd,char *line2);
113 enum CMDcode cmdpointsource(simptr sim,cmdptr cmd,char *line2);
114 enum CMDcode cmdvolumesource(simptr sim,cmdptr cmd,char *line2);
115 enum CMDcode cmdgaussiansource(simptr sim,cmdptr cmd,char *line2);
116 enum CMDcode cmdmovesurfacemol(simptr sim,cmdptr cmd,char *line2);
117 enum CMDcode cmdkillmol(simptr sim,cmdptr cmd,char *line2);
118 enum CMDcode cmdkillmolprob(simptr sim,cmdptr cmd,char *line2);
119 enum CMDcode cmdkillmolinsphere(simptr sim,cmdptr cmd,char *line2);
120 enum CMDcode cmdkillmolincmpt(simptr sim,cmdptr cmd,char *line2);
121 enum CMDcode cmdkillmoloutsidesystem(simptr sim,cmdptr cmd,char *line2);
122 enum CMDcode cmdfixmolcount(simptr sim,cmdptr cmd,char *line2);
123 enum CMDcode cmdfixmolcountrange(simptr sim,cmdptr cmd,char *line2);
124 enum CMDcode cmdfixmolcountonsurf(simptr sim,cmdptr cmd,char *line2);
125 enum CMDcode cmdfixmolcountrangeonsurf(simptr sim,cmdptr cmd,char *line2);
126 enum CMDcode cmdfixmolcountincmpt(simptr sim,cmdptr cmd,char *line2);
127 enum CMDcode cmdfixmolcountrangeincmpt(simptr sim,cmdptr cmd,char *line2);
128 enum CMDcode cmdequilmol(simptr sim,cmdptr cmd,char *line2);
129 enum CMDcode cmdreplacemol(simptr sim,cmdptr cmd,char *line2);
130 enum CMDcode cmdreplacexyzmol(simptr sim,cmdptr cmd,char *line2);
131 enum CMDcode cmdreplacevolmol(simptr sim,cmdptr cmd,char *line2);
132 enum CMDcode cmdreplacecmptmol(simptr sim,cmdptr cmd,char *line2);
133 enum CMDcode cmdmodulatemol(simptr sim,cmdptr cmd,char *line2);
134 enum CMDcode cmdreact1(simptr sim,cmdptr cmd,char *line2);
135 enum CMDcode cmdsetrateint(simptr sim,cmdptr cmd,char *line2);
136 enum CMDcode cmdshufflemollist(simptr sim,cmdptr cmd,char *line2);
137 enum CMDcode cmdshufflereactions(simptr sim,cmdptr cmd,char *line2);
138 enum CMDcode cmdsettimestep(simptr sim,cmdptr cmd,char *line2);
139 enum CMDcode cmdporttransport(simptr sim,cmdptr cmd,char *line2);
140 enum CMDcode cmdexcludebox(simptr sim,cmdptr cmd,char *line2);
141 enum CMDcode cmdexcludesphere(simptr sim,cmdptr cmd,char *line2);
142 enum CMDcode cmdincludeecoli(simptr sim,cmdptr cmd,char *line2);
143 enum CMDcode cmdsetreactionratemolcount(simptr sim,cmdptr cmd,char *line2);
144 enum CMDcode cmdexpandsystem(simptr sim,cmdptr cmd,char *line2);
145 enum CMDcode cmdtranslatecmpt(simptr sim,cmdptr cmd,char *line2);
146 enum CMDcode cmddiffusecmpt(simptr sim,cmdptr cmd,char *line2);
147 enum CMDcode cmdlongrangeforce(simptr sim,cmdptr cmd,char *line2);
148
149 // Smoldyn function declarations
150 double fnmolcount(void *voidsim,char *erstr,char *line2);
151 double fnmolcountonsurf(void *voidsim,char *erstr,char *line2);
152
153 // internal functions
154 void cmdv1free(cmdptr cmd);
155 void cmdv1v2free(cmdptr cmd);
156 void cmdListULVD4free(cmdptr cmd);
157 enum CMDcode conditionalcmdtype(simptr sim,cmdptr cmd,int nparam);
158 int insideecoli(double *pos,double *ofst,double rad,double length);
159 void putinecoli(double *pos,double *ofst,double rad,double length);
160 int molinpanels(simptr sim,moleculeptr mptr,int s,enum PanelShape ps);
161 void cmdmeansqrdispfree(cmdptr cmd);
162
163
164 /**********************************************************/
165 /********************* command processor ******************/
166 /**********************************************************/
167
168
169 /* docommand */
docommand(void * simvd,cmdptr cmd,char * line)170 enum CMDcode docommand(void *simvd,cmdptr cmd,char *line) {
171 simptr sim;
172 char word[STRCHAR],*line2;
173 int itct;
174
175 if(!simvd) return CMDok;
176 sim=(simptr) simvd;
177 if(!line) return CMDok;
178 itct=sscanf(line,"%s",word);
179 if(itct<=0) return CMDok;
180 line2=strnword(line,2);
181
182 Nvar=sim->nvar;
183 Varnames=sim->varnames;
184 Varvalues=sim->varvalues;
185
186 // simulation control
187 if(!strcmp(word,"stop")) return cmdstop(sim,cmd,line2);
188 else if(!strcmp(word,"pause")) return cmdpause(sim,cmd,line2);
189 else if(!strcmp(word,"beep")) return cmdbeep(sim,cmd,line2);
190 else if(!strcmp(word,"keypress")) return cmdkeypress(sim,cmd,line2);
191 else if(!strcmp(word,"setflag")) return cmdsetflag(sim,cmd,line2);
192 else if(!strcmp(word,"setrandseed")) return cmdsetrandseed(sim,cmd,line2);
193 else if(!strcmp(word,"setgraphics")) return cmdsetgraphics(sim,cmd,line2);
194 else if(!strcmp(word,"setgraphic_iter")) return cmdsetgraphic_iter(sim,cmd,line2);
195 else if(!strcmp(word,"updategraphics")) return cmdupdategraphics(sim,cmd,line2);
196
197 // file manipulation
198 else if(!strcmp(word,"overwrite")) return cmdoverwrite(sim,cmd,line2);
199 else if(!strcmp(word,"incrementfile")) return cmdincrementfile(sim,cmd,line2);
200
201 // conditional
202 else if(!strcmp(word,"ifflag")) return cmdifflag(sim,cmd,line2);
203 else if(!strcmp(word,"ifprob")) return cmdifprob(sim,cmd,line2);
204 else if(!strcmp(word,"ifno")) return cmdifno(sim,cmd,line2);
205 else if(!strcmp(word,"ifless")) return cmdifless(sim,cmd,line2);
206 else if(!strcmp(word,"ifmore")) return cmdifmore(sim,cmd,line2);
207 else if(!strcmp(word,"ifincmpt")) return cmdifincmpt(sim,cmd,line2);
208 else if(!strcmp(word,"ifchange")) return cmdifchange(sim,cmd,line2);
209 else if(!strcmp(word,"if")) return cmdif(sim,cmd,line2);
210
211 // system observation
212 else if(!strcmp(word,"echo")) return cmdecho(sim,cmd,line2);
213 else if(!strcmp(word,"evaluate")) return cmdevaluate(sim,cmd,line2);
214 else if(!strcmp(word,"warnescapee")) return cmdwarnescapee(sim,cmd,line2);
215 else if(!strcmp(word,"warnescapeecmpt")) return cmdwarnescapeecmpt(sim,cmd,line2);
216 else if(!strcmp(word,"molcountheader")) return cmdmolcountheader(sim,cmd,line2);
217 else if(!strcmp(word,"molcount")) return cmdmolcount(sim,cmd,line2);
218 else if(!strcmp(word,"molcountinbox")) return cmdmolcountinbox(sim,cmd,line2);
219 else if(!strcmp(word,"molcountincmpt")) return cmdmolcountincmpt(sim,cmd,line2);
220 else if(!strcmp(word,"molcountincmpts")) return cmdmolcountincmpts(sim,cmd,line2);
221 else if(!strcmp(word,"molcountincmpt2")) return cmdmolcountincmpt2(sim,cmd,line2);
222 else if(!strcmp(word,"molcountonsurf")) return cmdmolcountonsurf(sim,cmd,line2);
223 else if(!strcmp(word,"molcountspace")) return cmdmolcountspace(sim,cmd,line2);
224 else if(!strcmp(word,"molcountspace2d")) return cmdmolcountspace2d(sim,cmd,line2);
225 else if(!strcmp(word,"molcountspaceradial")) return cmdmolcountspaceradial(sim,cmd,line2);
226 else if(!strcmp(word,"molcountspacepolarangle")) return cmdmolcountspacepolarangle(sim,cmd,line2);
227 else if(!strcmp(word,"radialdistribution")) return cmdradialdistribution(sim,cmd,line2);
228 else if(!strcmp(word,"radialdistribution2")) return cmdradialdistribution2(sim,cmd,line2);
229 else if(!strcmp(word,"molcountspecies")) return cmdmolcountspecies(sim,cmd,line2);
230 else if(!strcmp(word,"molcountspecieslist")) return cmdmolcountspecieslist(sim,cmd,line2);
231 else if(!strcmp(word,"mollistsize")) return cmdmollistsize(sim,cmd,line2);
232 else if(!strcmp(word,"listmols")) return cmdlistmols(sim,cmd,line2);
233 else if(!strcmp(word,"listmols2")) return cmdlistmols2(sim,cmd,line2);
234 else if(!strcmp(word,"listmols3")) return cmdlistmols3(sim,cmd,line2);
235 else if(!strcmp(word,"listmols4")) return cmdlistmols4(sim,cmd,line2);
236 else if(!strcmp(word,"listmolscmpt")) return cmdlistmolscmpt(sim,cmd,line2);
237 else if(!strcmp(word,"molpos")) return cmdmolpos(sim,cmd,line2);
238 else if(!strcmp(word,"trackmol")) return cmdtrackmol(sim,cmd,line2);
239 else if(!strcmp(word,"molmoments")) return cmdmolmoments(sim,cmd,line2);
240 else if(!strcmp(word,"savesim")) return cmdsavesim(sim,cmd,line2);
241 else if(!strcmp(word,"meansqrdisp")) return cmdmeansqrdisp(sim,cmd,line2);
242 else if(!strcmp(word,"meansqrdisp2")) return cmdmeansqrdisp2(sim,cmd,line2);
243 else if(!strcmp(word,"meansqrdisp3")) return cmdmeansqrdisp3(sim,cmd,line2);
244 else if(!strcmp(word,"residencetime")) return cmdresidencetime(sim,cmd,line2);
245 else if(!strcmp(word,"diagnostics")) return cmddiagnostics(sim,cmd,line2);
246 else if(!strcmp(word,"executiontime")) return cmdexecutiontime(sim,cmd,line2);
247 else if(!strcmp(word,"writeVTK")) return cmdwriteVTK(sim,cmd,line2);
248 else if(!strcmp(word,"printLattice")) return cmdprintLattice(sim,cmd,line2);
249 else if(!strcmp(word,"printdata")) return cmdprintdata(sim,cmd,line2);
250
251 // system manipulation
252 else if(!strcmp(word,"set")) return cmdset(sim,cmd,line2);
253 else if(!strcmp(word,"pointsource")) return cmdpointsource(sim,cmd,line2);
254 else if(!strcmp(word,"volumesource")) return cmdvolumesource(sim,cmd,line2);
255 else if(!strcmp(word,"gaussiansource")) return cmdgaussiansource(sim,cmd,line2);
256 else if(!strcmp(word,"movesurfacemol")) return cmdmovesurfacemol(sim,cmd,line2);
257 else if(!strcmp(word,"killmol")) return cmdkillmol(sim,cmd,line2);
258 else if(!strcmp(word,"killmolprob")) return cmdkillmolprob(sim,cmd,line2);
259 else if(!strcmp(word,"killmolinsphere")) return cmdkillmolinsphere(sim,cmd,line2);
260 else if(!strcmp(word,"killmolincmpt")) return cmdkillmolincmpt(sim,cmd,line2);
261 else if(!strcmp(word,"killmoloutsidesystem")) return cmdkillmoloutsidesystem(sim,cmd,line2);
262 else if(!strcmp(word,"fixmolcount")) return cmdfixmolcount(sim,cmd,line2);
263 else if(!strcmp(word,"fixmolcountrange")) return cmdfixmolcountrange(sim,cmd,line2);
264 else if(!strcmp(word,"fixmolcountonsurf")) return cmdfixmolcountonsurf(sim,cmd,line2);
265 else if(!strcmp(word,"fixmolcountrangeonsurf")) return cmdfixmolcountrangeonsurf(sim,cmd,line2);
266 else if(!strcmp(word,"fixmolcountincmpt")) return cmdfixmolcountincmpt(sim,cmd,line2);
267 else if(!strcmp(word,"fixmolcountrangeincmpt")) return cmdfixmolcountrangeincmpt(sim,cmd,line2);
268 else if(!strcmp(word,"equilmol")) return cmdequilmol(sim,cmd,line2);
269 else if(!strcmp(word,"replacemol")) return cmdreplacemol(sim,cmd,line2);
270 else if(!strcmp(word,"replacexyzmol")) return cmdreplacexyzmol(sim,cmd,line2);
271 else if(!strcmp(word,"replacevolmol")) return cmdreplacevolmol(sim,cmd,line2);
272 else if(!strcmp(word,"replacecmptmol")) return cmdreplacecmptmol(sim,cmd,line2);
273 else if(!strcmp(word,"modulatemol")) return cmdmodulatemol(sim,cmd,line2);
274 else if(!strcmp(word,"react1")) return cmdreact1(sim,cmd,line2);
275 else if(!strcmp(word,"setrateint")) return cmdsetrateint(sim,cmd,line2);
276 else if(!strcmp(word,"shufflemollist")) return cmdshufflemollist(sim,cmd,line2);
277 else if(!strcmp(word,"shufflereactions")) return cmdshufflereactions(sim,cmd,line2);
278 else if(!strcmp(word,"settimestep")) return cmdsettimestep(sim,cmd,line2);
279 else if(!strcmp(word,"porttransport")) return cmdporttransport(sim,cmd,line2);
280 else if(!strcmp(word,"excludebox")) return cmdexcludebox(sim,cmd,line2);
281 else if(!strcmp(word,"excludesphere")) return cmdexcludesphere(sim,cmd,line2);
282 else if(!strcmp(word,"includeecoli")) return cmdincludeecoli(sim,cmd,line2);
283 else if(!strcmp(word,"setreactionratemolcount")) return cmdsetreactionratemolcount(sim,cmd,line2);
284 else if(!strcmp(word,"expandsystem")) return cmdexpandsystem(sim,cmd,line2);
285 else if(!strcmp(word,"translatecmpt")) return cmdtranslatecmpt(sim,cmd,line2);
286 else if(!strcmp(word,"diffusecmpt")) return cmddiffusecmpt(sim,cmd,line2);
287 else if(!strcmp(word,"longrangeforce")) return cmdlongrangeforce(sim,cmd,line2);
288
289 #ifdef VCELL
290 // vcell commands
291 else if(!strcmp(word,"vcellPrintProgress")) return cmdVCellPrintProgress(sim,cmd,line2);
292 else if(!strcmp(word,"vcellWriteOutput")) return cmdVCellWriteOutput(sim,cmd,line2);
293 else if(!strcmp(word,"vcellDataProcess")) return cmdVCellDataProcess(sim,cmd,line2);
294 #endif
295
296 SCMDCHECK(0,"command not recognized");
297 return CMDwarn; }
298
299
loadsmolfunctions(simptr sim)300 int loadsmolfunctions(simptr sim) {
301 double er;
302 char str1[STRCHAR],str2[STRCHAR];
303
304 er=0;
305 er+=strevalfunction(strcpy(str1,"molcount"),strcpy(str2,"dves"),(void*) sim,(void*) &fnmolcount,NULL,NULL,0);
306 er+=strevalfunction(strcpy(str1,"molcountonsurf"),strcpy(str2,"dves"),(void*) sim,(void*) &fnmolcountonsurf,NULL,NULL,0);
307
308 return (int) er; }
309
310
311 /**********************************************************/
312 /****************** simulation control ********************/
313 /**********************************************************/
314
315
316 /* cmdstop */
cmdstop(simptr sim,cmdptr cmd,char * line2)317 enum CMDcode cmdstop(simptr sim,cmdptr cmd,char *line2) {
318 if(line2 && !strcmp(line2,"cmdtype")) return CMDcontrol;
319 (void)sim;
320 (void)cmd;
321 return CMDstop; }
322
323
324 /* cmdpause */
cmdpause(simptr sim,cmdptr cmd,char * line2)325 enum CMDcode cmdpause(simptr sim,cmdptr cmd,char *line2) {
326 char c;
327 int tflag;
328
329 if(line2 && !strcmp(line2,"cmdtype")) return CMDcontrol;
330 if(!sim->graphss || sim->graphss->graphics==0) {
331 fprintf(stderr,"Simulation paused at time %g. Press enter to continue.",sim->time);
332 scanf("%c",&c);
333 return CMDok; }
334 tflag=strchr(sim->flags,'t')?1:0;
335 SCMDCHECK(sim->graphss && sim->graphss->graphics!=0 && !tflag,"pause doesn't work without graphics");
336 gl2State(1);
337 return CMDpause; }
338
339
340 /* cmdbeep */
cmdbeep(simptr sim,cmdptr cmd,char * line2)341 enum CMDcode cmdbeep(simptr sim,cmdptr cmd,char *line2) {
342 if(line2 && !strcmp(line2,"cmdtype")) return CMDcontrol;
343 (void)sim;
344 (void)cmd;
345 fprintf(stderr,"\7");
346 return CMDok; }
347
348
349 /* cmdkeypress */
cmdkeypress(simptr sim,cmdptr cmd,char * line2)350 enum CMDcode cmdkeypress(simptr sim,cmdptr cmd,char *line2) {
351 char c;
352 int itct,tflag;
353
354 if(line2 && !strcmp(line2,"cmdtype")) return CMDcontrol;
355 SCMDCHECK(line2,"missing argument");
356 itct=sscanf(line2,"%c",&c);
357 SCMDCHECK(itct==1,"cannot read character");
358 tflag=strchr(sim->flags,'t')?1:0;
359 SCMDCHECK(sim->graphss && sim->graphss->graphics!=0 && !tflag,"keypress doesn't work without graphics");
360 gl2SetKeyPush((unsigned char) c);
361 return CMDok; }
362
363
364 /* cmdsetflag */
cmdsetflag(simptr sim,cmdptr cmd,char * line2)365 enum CMDcode cmdsetflag(simptr sim,cmdptr cmd,char *line2) {
366 double f1;
367 int itct;
368
369 if(line2 && !strcmp(line2,"cmdtype")) return CMDcontrol;
370 SCMDCHECK(line2,"missing argument");
371 itct=strmathsscanf(line2,"%mlg",Varnames,Varvalues,Nvar,&f1);
372 SCMDCHECK(itct==1,"cannot read flag value");
373 scmdsetflag(sim->cmds,f1);
374 return CMDok; }
375
376
377 /* cmdsetrandseed */
cmdsetrandseed(simptr sim,cmdptr cmd,char * line2)378 enum CMDcode cmdsetrandseed(simptr sim,cmdptr cmd,char *line2) {
379 int itct;
380 long int seed;
381
382 if(line2 && !strcmp(line2,"cmdtype")) return CMDcontrol;
383 (void)sim;
384 SCMDCHECK(line2,"missing argument");
385 itct=sscanf(line2,"%li",&seed);
386 SCMDCHECK(itct==1,"cannot read seed");
387 if(seed<0) randomize((long int) time(NULL));
388 else randomize((long int) seed);
389 return CMDok; }
390
391
392 /* cmdsetgraphics */
cmdsetgraphics(simptr sim,cmdptr cmd,char * line2)393 enum CMDcode cmdsetgraphics(simptr sim,cmdptr cmd,char *line2) {
394 int itct;
395 char str[STRCHAR];
396
397 if(line2 && !strcmp(line2,"cmdtype")) return CMDcontrol;
398 if(!sim->graphss || sim->graphss->graphics==0) return CMDok;
399 SCMDCHECK(line2,"missing argument");
400 itct=sscanf(line2,"%s",str);
401 SCMDCHECK(itct==1,"cannot read graphics type");
402 if(!strcmp(str,"opengl")) sim->graphss->graphics=1;
403 else if(!strcmp(str,"opengl_good")) sim->graphss->graphics=2;
404 else SCMDCHECK(0,"unrecognized graphics type");
405 return CMDok; }
406
407 /* cmdsetgraphic_iter */
cmdsetgraphic_iter(simptr sim,cmdptr cmd,char * line2)408 enum CMDcode cmdsetgraphic_iter(simptr sim,cmdptr cmd,char *line2) {
409 int itct,iter;
410
411 if(line2 && !strcmp(line2,"cmdtype")) return CMDcontrol;
412 if(!sim->graphss || sim->graphss->graphics==0) return CMDok;
413 SCMDCHECK(line2,"missing argument");
414 itct=strmathsscanf(line2,"%mi",Varnames,Varvalues,Nvar,&iter);
415 SCMDCHECK(itct==1,"cannot read graphics iterations");
416 SCMDCHECK(iter>0,"graphics iterations must be >0");
417 sim->graphss->graphicit=iter;
418 return CMDok; }
419
420
421 /* cmdupdategraphics */
cmdupdategraphics(simptr sim,cmdptr cmd,char * line2)422 enum CMDcode cmdupdategraphics(simptr sim,cmdptr cmd,char *line2) {
423 if(line2 && !strcmp(line2,"cmdtype")) return CMDcontrol;
424 (void)cmd;
425 if(!sim->graphss || sim->graphss->graphics==0) return CMDok;
426 smolPostRedisplay();
427 return CMDok; }
428
429
430 /**********************************************************/
431 /****************** file manipulation ********************/
432 /**********************************************************/
433
434
435 /* cmdoverwrite */
cmdoverwrite(simptr sim,cmdptr cmd,char * line2)436 enum CMDcode cmdoverwrite(simptr sim,cmdptr cmd,char *line2) {
437 int er;
438
439 if(line2 && !strcmp(line2,"cmdtype")) return CMDcontrol;
440 SCMDCHECK(line2,"missing argument");
441 er=scmdoverwrite(sim->cmds,line2);
442 SCMDCHECK(er!=1,"file not declared");
443 SCMDCHECK(er!=2,"failed to open file for writing");
444 return CMDok; }
445
446
447 /* cmdincrementfile */
cmdincrementfile(simptr sim,cmdptr cmd,char * line2)448 enum CMDcode cmdincrementfile(simptr sim,cmdptr cmd,char *line2) {
449 int er;
450
451 if(line2 && !strcmp(line2,"cmdtype")) return CMDcontrol;
452 SCMDCHECK(line2,"missing argument");
453 er=scmdincfile(sim->cmds,line2);
454 SCMDCHECK(er!=1,"file name not declared");
455 SCMDCHECK(er!=2,"failed to open new file for writing");
456 return CMDok; }
457
458
459
460 /**********************************************************/
461 /********************** conditional ***********************/
462 /**********************************************************/
463
464
465 /* cmdifflag */
cmdifflag(simptr sim,cmdptr cmd,char * line2)466 enum CMDcode cmdifflag(simptr sim,cmdptr cmd,char *line2) {
467 int itct;
468 char ch;
469 double f1,flag;
470
471 if(line2 && !strcmp(line2,"cmdtype")) return conditionalcmdtype(sim,cmd,2);
472 SCMDCHECK(line2,"missing arguments");
473 itct=strmathsscanf(line2,"%c %mlg",Varnames,Varvalues,Nvar,&ch,&f1);
474 SCMDCHECK(itct==2,"cannot read comparison symbol or flag value");
475 SCMDCHECK(ch=='<' || ch=='=' || ch=='>',"comparison symbol has to be <, =, or >");
476 flag=scmdreadflag(sim->cmds);
477 if((ch=='<' && flag<f1) || (ch=='=' && flag==f1) || (ch=='>' && flag>f1))
478 return docommand(sim,cmd,strnword(line2,3));
479 return CMDok; }
480
481
482 /* cmdifprob */
cmdifprob(simptr sim,cmdptr cmd,char * line2)483 enum CMDcode cmdifprob(simptr sim,cmdptr cmd,char *line2) {
484 int itct;
485 double f1;
486
487 if(line2 && !strcmp(line2,"cmdtype")) return conditionalcmdtype(sim,cmd,1);
488 SCMDCHECK(line2,"missing arguments");
489 itct=strmathsscanf(line2,"%mlg",Varnames,Varvalues,Nvar,&f1);
490 SCMDCHECK(itct==1,"cannot read probability value");
491 SCMDCHECK(f1>=0 && f1<=1,"probability value needs to be between 0 and 1");
492 if(randCOD()<f1)
493 return docommand(sim,cmd,strnword(line2,2));
494 return CMDok; }
495
496
497 /* cmdifno */
cmdifno(simptr sim,cmdptr cmd,char * line2)498 enum CMDcode cmdifno(simptr sim,cmdptr cmd,char *line2) {
499 int i,count,*index;
500 enum MolecState ms;
501
502 if(line2 && !strcmp(line2,"cmdtype")) {
503 return conditionalcmdtype(sim,cmd,1); }
504
505 i=molstring2index1(sim,line2,&ms,&index);
506 SCMDCHECK(i!=-1,"species is missing or cannot be read");
507 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
508 SCMDCHECK(i!=-3,"cannot read molecule state value");
509 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
510 SCMDCHECK(i!=-7,"error allocating memory");
511 count=(i==-4)?0:molcount(sim,i,index,ms,1);
512 if(count) return CMDok;
513 return docommand(sim,cmd,strnword(line2,2)); }
514
515
516 /* cmdifless */
cmdifless(simptr sim,cmdptr cmd,char * line2)517 enum CMDcode cmdifless(simptr sim,cmdptr cmd,char *line2) {
518 int itct,i,*index,count,min;
519 enum MolecState ms;
520
521 if(line2 && !strcmp(line2,"cmdtype")) {
522 return conditionalcmdtype(sim,cmd,2); }
523
524 i=molstring2index1(sim,line2,&ms,&index);
525 SCMDCHECK(i!=-1,"species is missing or cannot be read");
526 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
527 SCMDCHECK(i!=-3,"cannot read molecule state value");
528 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
529 SCMDCHECK(i!=-7,"error allocating memory");
530 SCMDCHECK(line2=strnword(line2,2),"missing value argument");
531 itct=strmathsscanf(line2,"%mi",Varnames,Varvalues,Nvar,&min);
532 SCMDCHECK(itct==1,"cannot read value argument");
533 count=(i==-4)?0:molcount(sim,i,index,ms,min);
534 if(count<min) return docommand(sim,cmd,strnword(line2,2));
535 return CMDok; }
536
537
538 /* cmdifmore */
cmdifmore(simptr sim,cmdptr cmd,char * line2)539 enum CMDcode cmdifmore(simptr sim,cmdptr cmd,char *line2) {
540 int itct,i,*index,count,min;
541 enum MolecState ms;
542
543 if(line2 && !strcmp(line2,"cmdtype")) {
544 return conditionalcmdtype(sim,cmd,2); }
545
546 i=molstring2index1(sim,line2,&ms,&index);
547 SCMDCHECK(i!=-1,"species is missing or cannot be read");
548 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
549 SCMDCHECK(i!=-3,"cannot read molecule state value");
550 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
551 SCMDCHECK(i!=-7,"error allocating memory");
552 SCMDCHECK(line2=strnword(line2,2),"missing value argument");
553 itct=strmathsscanf(line2,"%mi",Varnames,Varvalues,Nvar,&min);
554 SCMDCHECK(itct==1,"cannot read value argument");
555 count=(i==-4)?0:molcount(sim,i,index,ms,min+1);
556 if(count>min) return docommand(sim,cmd,strnword(line2,2));
557 return CMDok; }
558
559
560 /* cmdifincmpt */
cmdifincmpt(simptr sim,cmdptr cmd,char * line2)561 enum CMDcode cmdifincmpt(simptr sim,cmdptr cmd,char *line2) {
562 int itct,i,min,c,*index;
563 enum MolecState ms;
564 char cname[STRCHAR];
565 compartssptr cmptss;
566 char ch;
567 moleculeptr mptr;
568 static compartptr cmpt=NULL;
569 static int inscan=0,count=0;
570
571 if(inscan) goto scanportion;
572 if(line2 && !strcmp(line2,"cmdtype")) return conditionalcmdtype(sim,cmd,4);
573
574 cmptss=sim->cmptss;
575 SCMDCHECK(cmptss,"no compartments defined");
576 SCMDCHECK(sim->mols,"molecules are undefined");
577 SCMDCHECK(line2,"missing argument");
578 i=molstring2index1(sim,line2,&ms,&index);
579 SCMDCHECK(i!=-1,"species is missing or cannot be read");
580 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
581 SCMDCHECK(i!=-3,"cannot read molecule state value");
582 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
583 SCMDCHECK(i!=-7,"error allocating memory");
584 SCMDCHECK(line2=strnword(line2,2),"missing value argument");
585 itct=strmathsscanf(line2,"%c %mi %s",Varnames,Varvalues,Nvar,&ch,&min,cname);
586 SCMDCHECK(itct==3,"cannot read symbol, value, and/or compartment arguments");
587 SCMDCHECK(ch=='<' || ch=='=' || ch=='>',"comparison symbol has to be <, =, or >");
588 c=stringfind(cmptss->cnames,cmptss->ncmpt,cname);
589 SCMDCHECK(c>=0,"compartment name not recognized");
590 cmpt=cmptss->cmptlist[c];
591 line2=strnword(line2,4);
592
593 if(i==-4) count=0;
594 else {
595 count=0;
596 inscan=1;
597 molscancmd(sim,i,index,ms,cmd,cmdifincmpt);
598 inscan=0; }
599 if((ch=='<' && count<min) || (ch=='=' && count==min) || (ch=='>' && count>min))
600 return docommand(sim,cmd,line2);
601 return CMDok;
602
603 scanportion:
604 mptr=(moleculeptr) line2;
605 if(posincompart(sim,mptr->pos,cmpt,0)) count++;
606 return CMDok; }
607
608
609 /* cmdifchange */
cmdifchange(simptr sim,cmdptr cmd,char * line2)610 enum CMDcode cmdifchange(simptr sim,cmdptr cmd,char *line2) {
611 int itct,i,*index,count,num,diff;
612 enum MolecState ms;
613 char change;
614
615 if(line2 && !strcmp(line2,"cmdtype")) {
616 return conditionalcmdtype(sim,cmd,2); }
617
618 i=molstring2index1(sim,line2,&ms,&index);
619 SCMDCHECK(i!=-1,"species is missing or cannot be read");
620 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
621 SCMDCHECK(i!=-3,"cannot read molecule state value");
622 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
623 SCMDCHECK(i!=-7,"error allocating memory");
624 SCMDCHECK(line2=strnword(line2,2),"missing value argument");
625 itct=strmathsscanf(line2,"%c %mi",Varnames,Varvalues,Nvar,&change,&num);
626 SCMDCHECK(itct==2,"cannot read change or number arguments");
627 SCMDCHECK(line2=strnword(line2,3),"missing conditional command");
628
629 if(cmd->i1==0) {
630 cmd->i1=1;
631 cmd->i2=(i==-4)?0:molcount(sim,i,index,ms,-1); }
632 else {
633 count=(i==-4)?0:molcount(sim,i,index,ms,-1);
634 diff=count-cmd->i2;
635 cmd->i2=count;
636 if((change=='>' && diff>num) || (change=='<' && diff<num) || (change=='=' && diff==num) || (change=='!' && diff!=num))
637 return docommand(sim,cmd,line2); }
638 return CMDok; }
639
640
641 /* cmdif */
cmdif(simptr sim,cmdptr cmd,char * line2)642 enum CMDcode cmdif(simptr sim,cmdptr cmd,char *line2) {
643 int itct;
644 char symbol;
645 double value1,value2;
646
647 if(line2 && !strcmp(line2,"cmdtype")) {
648 return conditionalcmdtype(sim,cmd,2); }
649
650 itct=strmathsscanf(line2,"%mlg %c %mlg",Varnames,Varvalues,Nvar,&value1,&symbol,&value2);
651 SCMDCHECK(itct==3,"cannot read command arguments");
652 SCMDCHECK(line2=strnword(line2,4),"missing conditional command");
653
654 if((symbol=='>' && value1>value2) || (symbol=='<' && value1<value2) || (symbol=='=' && value1==value2))
655 return docommand(sim,cmd,line2);
656
657 return CMDok; }
658
659
660 /**********************************************************/
661 /***************** observation commands *******************/
662 /**********************************************************/
663
664
665 /* cmdwarnescapee */
cmdwarnescapee(simptr sim,cmdptr cmd,char * line2)666 enum CMDcode cmdwarnescapee(simptr sim,cmdptr cmd,char *line2) {
667 int i,escape,*index,er;
668 enum MolecState ms;
669 moleculeptr mptr;
670 double *pos,*posx,*via;
671 char string[STRCHAR];
672 static int inscan=0;
673 static FILE *fptr=NULL;
674 static int dataid=-1;
675
676 if(inscan) goto scanportion;
677 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
678
679 i=molstring2index1(sim,line2,&ms,&index);
680 SCMDCHECK(i!=-1,"species is missing or cannot be read");
681 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
682 SCMDCHECK(i!=-3,"cannot read molecule state value");
683 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
684 SCMDCHECK(i!=-7,"error allocating memory");
685 line2=strnword(line2,2);
686 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
687 SCMDCHECK(er!=-1,"file or data name not recognized");
688
689 if(i!=-4) {
690 inscan=1;
691 molscancmd(sim,i,index,ms,cmd,cmdwarnescapee);
692 inscan=0;
693 scmdflush(fptr); }
694 return CMDok;
695
696 scanportion:
697 mptr=(moleculeptr)line2;
698 pos=mptr->pos;
699 escape=!posinsystem(sim,pos);
700 if(escape) {
701 posx=mptr->posx;
702 escape=!posinsystem(sim,posx);
703 if(!escape) {
704 via=mptr->via;
705 if(sim->dim==1) {
706 scmdfprintf(cmd->cmds,fptr,"New escapee: %g #%s %g to %g via %g\n",sim->time,molserno2string(mptr->serno,string),posx[0],pos[0],via[0]);
707 scmdappenddata(cmd->cmds,dataid,1,5,sim->time,(double)(mptr->serno),posx[0],pos[0],via[0]); }
708 else if(sim->dim==2) {
709 scmdfprintf(cmd->cmds,fptr,"New escapee: %g #%s (%g,%g) to (%g,%g) via (%g,%g)\n",sim->time,molserno2string(mptr->serno,string),posx[0],posx[1],pos[0],pos[1],via[0],via[1]);
710 scmdappenddata(cmd->cmds,dataid,1,8,sim->time,(double)(mptr->serno),posx[0],posx[1],pos[0],pos[1],via[0],via[1]); }
711 else {
712 scmdfprintf(cmd->cmds,fptr,"New escapee: %g #%s (%g,%g,%g) to (%g,%g,%g) via (%g,%g,%g)\n",sim->time,molserno2string(mptr->serno,string),posx[0],posx[1],posx[2],pos[0],pos[1],pos[2],via[0],via[1],via[2]);
713 scmdappenddata(cmd->cmds,dataid,1,11,sim->time,(double)(mptr->serno),posx[0],posx[1],posx[2],pos[0],pos[1],pos[2],via[0],via[1],via[2]); }}}
714 return CMDok; }
715
716
717 /* cmdwarnescapeecmpt */
cmdwarnescapeecmpt(simptr sim,cmdptr cmd,char * line2)718 enum CMDcode cmdwarnescapeecmpt(simptr sim,cmdptr cmd,char *line2) {
719 int i,escape,*index,er,itct,c;
720 enum MolecState ms;
721 moleculeptr mptr;
722 compartssptr cmptss;
723 double *pos,*posx,*via;
724 char string[STRCHAR],nm[STRCHAR];
725 static int inscan=0;
726 static compartptr cmpt=NULL;
727 static FILE *fptr=NULL;
728 static int dataid=-1;
729
730 if(inscan) goto scanportion;
731 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
732
733 i=molstring2index1(sim,line2,&ms,&index);
734 SCMDCHECK(i!=-1,"species is missing or cannot be read");
735 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
736 SCMDCHECK(i!=-3,"cannot read molecule state value");
737 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
738 SCMDCHECK(i!=-7,"error allocating memory");
739 line2=strnword(line2,2);
740 cmptss=sim->cmptss;
741 SCMDCHECK(cmptss,"no compartments defined");
742 SCMDCHECK(line2,"missing argument");
743 itct=sscanf(line2,"%s",nm);
744 SCMDCHECK(itct==1,"cannot read argument");
745 c=stringfind(cmptss->cnames,cmptss->ncmpt,nm);
746 SCMDCHECK(c>=0,"compartment name not recognized");
747 cmpt=cmptss->cmptlist[c];
748 line2=strnword(line2,2);
749 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
750 SCMDCHECK(er!=-1,"file or data name not recognized");
751
752 if(i!=-4) {
753 inscan=1;
754 molscancmd(sim,i,index,ms,cmd,cmdwarnescapeecmpt);
755 inscan=0;
756 scmdflush(fptr); }
757 return CMDok;
758
759 scanportion:
760 mptr=(moleculeptr)line2;
761 pos=mptr->pos;
762 escape=!posincompart(sim,pos,cmpt,0);
763 if(escape) {
764 posx=mptr->posx;
765 escape=!posincompart(sim,posx,cmpt,0);
766 if(!escape) {
767 via=mptr->via;
768 if(sim->dim==1) {
769 scmdfprintf(cmd->cmds,fptr,"New escapee: %g #%s %g to %g via %g\n",sim->time,molserno2string(mptr->serno,string),posx[0],pos[0],via[0]);
770 scmdappenddata(cmd->cmds,dataid,1,5,sim->time,(double)(mptr->serno),posx[0],pos[0],via[0]); }
771 else if(sim->dim==2) {
772 scmdfprintf(cmd->cmds,fptr,"New escapee: %g #%s (%g,%g) to (%g,%g) via (%g,%g)\n",sim->time,molserno2string(mptr->serno,string),posx[0],posx[1],pos[0],pos[1],via[0],via[1]);
773 scmdappenddata(cmd->cmds,dataid,1,8,sim->time,(double)(mptr->serno),posx[0],posx[1],pos[0],pos[1],via[0],via[1]); }
774 else {
775 scmdfprintf(cmd->cmds,fptr,"New escapee: %g #%s (%g,%g,%g) to (%g,%g,%g) via (%g,%g,%g)\n",sim->time,molserno2string(mptr->serno,string),posx[0],posx[1],posx[2],pos[0],pos[1],pos[2],via[0],via[1],via[2]);
776 scmdappenddata(cmd->cmds,dataid,1,11,sim->time,(double)(mptr->serno),posx[0],posx[1],posx[2],pos[0],pos[1],pos[2],via[0],via[1],via[2]); }}}
777 return CMDok; }
778
779
780 /* cmdecho */
cmdecho(simptr sim,cmdptr cmd,char * line2)781 enum CMDcode cmdecho(simptr sim,cmdptr cmd,char *line2) {
782 int er;
783 FILE *fptr;
784 char *termqt,str[STRCHAR];
785
786 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
787 er=scmdgetfptr(sim->cmds,line2,1,&fptr,NULL);
788 SCMDCHECK(er!=-1,"file name not recognized");
789 line2=strnword(line2,2);
790 SCMDCHECK(line2=strchr(line2,'"'),"no starting quote on string");
791 strncpy(str,line2+1,STRCHAR-1);
792 str[STRCHAR-1]='\0';
793 SCMDCHECK(termqt=strchr(str,'"'),"no terminal quote on string");
794 *termqt='\0';
795 strbslash2escseq(str);
796 scmdfprintf(cmd->cmds,fptr,"%s",str);
797 scmdflush(fptr);
798 return CMDok; }
799
800
801 /* cmdevaluate */
cmdevaluate(simptr sim,cmdptr cmd,char * line2)802 enum CMDcode cmdevaluate(simptr sim,cmdptr cmd,char *line2) {
803 FILE *fptr;
804 double answer;
805 int it,er,dataid;
806 char erstr[STRCHAR];
807
808 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
809 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
810 SCMDCHECK(er!=-1,"file or data name not recognized");
811 line2=strnword(line2,2);
812 SCMDCHECK(line2,"missing item to evaluate");
813 it=strmathsscanf(line2,"%mlg",sim->varnames,sim->varvalues,sim->nvar,&answer);
814 if(it!=1) {
815 er=strmatherror(erstr,1);
816 SCMDCHECK(!er,"%s",erstr); }
817 scmdfprintf(cmd->cmds,fptr,"%g\n",answer);
818 scmdappenddata(cmd->cmds,dataid,1,1,answer);
819 scmdflush(fptr);
820 return CMDok; }
821
822
823 /* cmdmolcountheader */
cmdmolcountheader(simptr sim,cmdptr cmd,char * line2)824 enum CMDcode cmdmolcountheader(simptr sim,cmdptr cmd,char *line2) {
825 FILE *fptr;
826 int i,er;
827
828 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
829 er=scmdgetfptr(sim->cmds,line2,1,&fptr,NULL);
830 SCMDCHECK(er!=-1,"file name not recognized");
831 SCMDCHECK(sim->mols,"molecules are undefined");
832 scmdfprintf(cmd->cmds,fptr,"time");
833 for(i=1;i<sim->mols->nspecies;i++)
834 scmdfprintf(cmd->cmds,fptr,"%,%s",sim->mols->spname[i]);
835 scmdfprintf(cmd->cmds,fptr,"\n");
836 scmdflush(fptr);
837 return CMDok; }
838
839
840 /* cmdmolcount */
cmdmolcount(simptr sim,cmdptr cmd,char * line2)841 enum CMDcode cmdmolcount(simptr sim,cmdptr cmd,char *line2) {
842 FILE *fptr;
843 int i,nspecies,*ctlat,ilat,er,dataid;
844 latticeptr lat;
845 moleculeptr mptr;
846 static int inscan=0,*ct=NULL;
847
848 if(inscan) goto scanportion;
849 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
850
851 SCMDCHECK(cmd->i1!=-1,"error on setup"); // failed before, don't try again
852 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
853 SCMDCHECK(er!=-1,"file or data name not recognized");
854 SCMDCHECK(sim->mols,"molecules are undefined");
855
856 nspecies=sim->mols->nspecies;
857 if(cmd->i1!=nspecies) { // allocate counter if required
858 cmdv1free(cmd);
859 cmd->i1=nspecies;
860 cmd->freefn=&cmdv1v2free;
861 cmd->v1=calloc(nspecies,sizeof(int));
862 if(!cmd->v1) {cmd->i1=-1;return CMDwarn;} }
863
864 ct=(int*)cmd->v1;
865 for(i=0;i<nspecies;i++) ct[i]=0; // clear counters
866
867 inscan=1;
868 molscancmd(sim,-1,NULL,MSall,cmd,cmdmolcount);
869 inscan=0;
870
871 if(sim->latticess) {
872 if(cmd->i2!=nspecies) {
873 free(cmd->v2);
874 cmd->i2=nspecies;
875 cmd->v2=calloc(nspecies,sizeof(int));
876 if(!cmd->v2) {cmd->i2=-1;return CMDwarn;}}
877
878 ctlat=(int*)cmd->v2;
879 for(ilat=0;ilat<sim->latticess->nlattice;++ilat) {
880 lat=sim->latticess->latticelist[ilat];
881 for(i=1;i<nspecies;i++) ctlat[i]=0;
882 if(lat->type==LATTICEnsv) {
883 NSV_CALL(nsv_molcount(lat->nsv,ctlat)); }
884 else if(lat->type==LATTICEpde) {
885 //not implemented
886 }
887 for(i=1;i<nspecies;i++) {
888 ct[i]+=ctlat[i]; }}}
889
890 scmdfprintf(cmd->cmds,fptr,"%g",sim->time);
891 scmdappenddata(cmd->cmds,dataid,1,1,sim->time);
892 for(i=1;i<nspecies;i++) {
893 scmdfprintf(cmd->cmds,fptr,"%,%i",ct[i]);
894 scmdappenddata(cmd->cmds,dataid,0,1,(double)ct[i]); }
895 scmdfprintf(cmd->cmds,fptr,"\n");
896 scmdflush(fptr);
897 return CMDok;
898
899 scanportion:
900 mptr=(moleculeptr) line2;
901 ct[mptr->ident]++;
902 return CMDok; }
903
904
905 /* fnmolcount */
fnmolcount(void * voidsim,char * erstr,char * line2)906 double fnmolcount(void *voidsim,char *erstr,char *line2) {
907 simptr sim;
908 enum MolecState ms;
909 int i,*index;
910 static char oldline2[STRCHAR]="\0";
911 static int inscan=0,ct;
912 static long int oldtouch=0;
913
914 sim=(simptr) voidsim;
915 if(inscan) goto scanportion;
916
917 if(!sim->mols) return 0;
918 if(sim->mols->touch==oldtouch && !strcmp(line2,oldline2)) return ct;
919 strcpy(oldline2,line2);
920 oldtouch=sim->mols->touch;
921
922 SFNCHECK(line2,"missing arguments");
923 i=molstring2index1(sim,line2,&ms,&index);
924 SFNCHECK(i!=-1,"species is missing or cannot be read");
925 SFNCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
926 SFNCHECK(i!=-3,"cannot read molecule state value");
927 SFNCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
928 SFNCHECK(i!=-7,"error allocating memory");
929
930 ct=0;
931 inscan=1;
932 molscanfn(sim,i,index,ms,erstr,fnmolcount);
933 inscan=0;
934 return ct;
935
936 scanportion:
937 ct++;
938 return 0; }
939
940
941 /* fnmolcountonsurf */
fnmolcountonsurf(void * voidsim,char * erstr,char * line2)942 double fnmolcountonsurf(void *voidsim,char *erstr,char *line2) {
943 simptr sim;
944 enum MolecState ms;
945 int i,*index,s,comma,itct;
946 static int inscan=0,ct;
947 surfacessptr srfss;
948 char nm[STRCHAR];
949 static surfaceptr srf;
950 moleculeptr mptr;
951 static long int oldtouch=0;
952 static char oldline2[STRCHAR]="\0";
953
954 sim=(simptr) voidsim;
955 if(inscan) goto scanportion;
956
957 if(!sim->mols) return 0;
958 if(sim->mols->touch==oldtouch && !strcmp(line2,oldline2)) return ct;
959 strcpy(oldline2,line2);
960 oldtouch=sim->mols->touch;
961
962 srfss=sim->srfss;
963 SFNCHECK(srfss,"no surfaces defined");
964 SFNCHECK(line2,"missing arguments");
965 SFNCHECK((comma=strChrBrackets(line2,-1,',',"([{,\"'"))>0,"missing parameter");
966 line2[comma]='\0';
967 i=molstring2index1(sim,line2,&ms,&index);
968 SFNCHECK(i!=-1,"species is missing or cannot be read");
969 SFNCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
970 SFNCHECK(i!=-3,"cannot read molecule state value");
971 SFNCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
972 SFNCHECK(i!=-7,"error allocating memory");
973 line2+=comma+1;
974 itct=sscanf(line2,"%s",nm);
975 SFNCHECK(itct==1,"cannot read surface name");
976 s=stringfind(srfss->snames,srfss->nsrf,nm);
977 SFNCHECK(s>=0,"surface name '%s' not recognized",nm);
978 srf=srfss->srflist[s];
979
980 ct=0;
981 inscan=1;
982 molscanfn(sim,i,index,ms,erstr,fnmolcountonsurf);
983 inscan=0;
984 return ct;
985
986 scanportion: //?? This is very inefficient; should use surface molecule list
987 mptr=(moleculeptr) line2;
988 if(mptr->mstate!=MSsoln && mptr->pnl->srf==srf) ct++;
989 return 0; }
990
991
992 /* cmdmolcountinbox */
cmdmolcountinbox(simptr sim,cmdptr cmd,char * line2)993 enum CMDcode cmdmolcountinbox(simptr sim,cmdptr cmd,char *line2) {
994 FILE *fptr;
995 int d,dim,itct,i,nspecies,er,dataid;
996 moleculeptr mptr;
997 static double low[3]={0,0,0},high[3]={0,0,0};
998 static int inscan=0,*ct;
999
1000 if(inscan) goto scanportion;
1001 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
1002
1003 SCMDCHECK(cmd->i1!=-1,"error on setup"); // failed before, don't try again
1004 SCMDCHECK(sim->mols,"molecules are undefined");
1005 dim=sim->dim;
1006 for(d=0;d<dim;d++) {
1007 SCMDCHECK(line2,"missing argument");
1008 itct=strmathsscanf(line2,"%mlg %mlg",Varnames,Varvalues,Nvar,&low[d],&high[d]);
1009 SCMDCHECK(itct==2,"read failure");
1010 line2=strnword(line2,3); }
1011 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
1012 SCMDCHECK(er!=-1,"file or data name not recognized");
1013
1014 nspecies=sim->mols->nspecies;
1015 if(cmd->i1!=nspecies) { // allocate counter if required
1016 cmdv1free(cmd);
1017 cmd->i1=nspecies;
1018 cmd->freefn=&cmdv1free;
1019 cmd->v1=calloc(nspecies,sizeof(int));
1020 if(!cmd->v1) {cmd->i1=-1;return CMDwarn;} }
1021
1022 ct=(int*)cmd->v1;
1023 for(i=0;i<nspecies;i++) ct[i]=0;
1024
1025 inscan=1;
1026 molscancmd(sim,-1,NULL,MSall,cmd,cmdmolcountinbox);
1027 inscan=0;
1028
1029 scmdfprintf(cmd->cmds,fptr,"%g",sim->time);
1030 scmdappenddata(cmd->cmds,dataid,1,1,sim->time);
1031 for(i=1;i<nspecies;i++) {
1032 scmdfprintf(cmd->cmds,fptr,"%,%i",ct[i]);
1033 scmdappenddata(cmd->cmds,dataid,0,1,(double)ct[i]); }
1034 scmdfprintf(cmd->cmds,fptr,"\n");
1035 scmdflush(fptr);
1036 return CMDok;
1037
1038 scanportion:
1039 mptr=(moleculeptr) line2;
1040 for(d=0;d<sim->dim;d++)
1041 if(mptr->pos[d]<low[d] || mptr->pos[d]>high[d]) return CMDok;
1042 ct[mptr->ident]++;
1043 return CMDok; }
1044
1045
1046 /* cmdmolcountincmpt */
cmdmolcountincmpt(simptr sim,cmdptr cmd,char * line2)1047 enum CMDcode cmdmolcountincmpt(simptr sim,cmdptr cmd,char *line2) {
1048 FILE *fptr;
1049 char nm[STRCHAR];
1050 compartssptr cmptss;
1051 int itct,c,i,nspecies,er,dataid;
1052 moleculeptr mptr;
1053 static compartptr cmpt=NULL;
1054 static int inscan=0,*ct;
1055
1056 if(inscan) goto scanportion;
1057 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
1058
1059 SCMDCHECK(cmd->i1!=-1,"error on setup"); // failed before, don't try again
1060 cmptss=sim->cmptss;
1061 SCMDCHECK(cmptss,"no compartments defined");
1062 SCMDCHECK(sim->mols,"molecules are undefined");
1063 SCMDCHECK(line2,"missing argument");
1064 itct=sscanf(line2,"%s",nm);
1065 SCMDCHECK(itct==1,"cannot read argument");
1066 c=stringfind(cmptss->cnames,cmptss->ncmpt,nm);
1067 SCMDCHECK(c>=0,"compartment name not recognized");
1068 cmpt=cmptss->cmptlist[c];
1069 line2=strnword(line2,2);
1070 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
1071 SCMDCHECK(er!=-1,"file or data name not recognized");
1072
1073 nspecies=sim->mols->nspecies;
1074 if(cmd->i1!=nspecies) { // allocate counter if required
1075 cmdv1free(cmd);
1076 cmd->i1=nspecies;
1077 cmd->freefn=&cmdv1free;
1078 cmd->v1=calloc(nspecies,sizeof(int));
1079 if(!cmd->v1) {cmd->i1=-1;return CMDwarn;} }
1080
1081 ct=(int*)cmd->v1;
1082 for(i=0;i<nspecies;i++) ct[i]=0;
1083
1084 inscan=1;
1085 molscancmd(sim,-1,NULL,MSsoln,cmd,cmdmolcountincmpt);
1086 inscan=0;
1087
1088 scmdfprintf(cmd->cmds,fptr,"%g",sim->time);
1089 scmdappenddata(cmd->cmds,dataid,1,1,sim->time);
1090 for(i=1;i<nspecies;i++) {
1091 scmdfprintf(cmd->cmds,fptr,"%,%i",ct[i]);
1092 scmdappenddata(cmd->cmds,dataid,0,1,(double)ct[i]); }
1093 scmdfprintf(cmd->cmds,fptr,"\n");
1094 scmdflush(fptr);
1095 return CMDok;
1096
1097 scanportion:
1098 mptr=(moleculeptr) line2;
1099 if(posincompart(sim,mptr->pos,cmpt,0)) ct[mptr->ident]++;
1100 return CMDok; }
1101
1102
1103 /* cmdmolcountincmpts */
cmdmolcountincmpts(simptr sim,cmdptr cmd,char * line2)1104 enum CMDcode cmdmolcountincmpts(simptr sim,cmdptr cmd,char *line2) {
1105 FILE *fptr;
1106 char nm[STRCHAR];
1107 compartssptr cmptss;
1108 int itct,c,i,ic,er,dataid;
1109 moleculeptr mptr;
1110 static int inscan=0,*ct,ncmpt,nspecies;
1111 static compartptr cmptlist[16];
1112
1113 if(inscan) goto scanportion;
1114 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
1115
1116 SCMDCHECK(cmd->i1!=-1,"error on setup"); // failed before, don't try again
1117 cmptss=sim->cmptss;
1118 SCMDCHECK(cmptss,"no compartments defined");
1119 SCMDCHECK(sim->mols,"molecules are undefined");
1120 SCMDCHECK(line2,"missing argument");
1121 ncmpt=wordcount(line2)-1;
1122 SCMDCHECK(ncmpt>=1,"no compartment or no output file listed");
1123 for(ic=0;ic<ncmpt;ic++) {
1124 itct=sscanf(line2,"%s",nm);
1125 SCMDCHECK(itct==1,"cannot read compartment name");
1126 c=stringfind(cmptss->cnames,cmptss->ncmpt,nm);
1127 SCMDCHECK(c>=0,"compartment name not recognized");
1128 cmptlist[ic]=cmptss->cmptlist[c];
1129 line2=strnword(line2,2);
1130 SCMDCHECK(line2,"missing argument"); }
1131 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
1132 SCMDCHECK(er!=-1,"file or data name not recognized");
1133
1134 nspecies=sim->mols->nspecies;
1135 if(cmd->i1!=nspecies) { // allocate counter if required
1136 cmdv1free(cmd);
1137 cmd->i1=nspecies;
1138 cmd->freefn=&cmdv1free;
1139 cmd->v1=calloc(nspecies*ncmpt,sizeof(int));
1140 if(!cmd->v1) {cmd->i1=-1;return CMDwarn;} }
1141
1142 ct=(int*)cmd->v1;
1143 for(i=0;i<nspecies*ncmpt;i++) ct[i]=0;
1144
1145 inscan=1;
1146 molscancmd(sim,-1,NULL,MSsoln,cmd,cmdmolcountincmpts);
1147 inscan=0;
1148
1149 scmdfprintf(cmd->cmds,fptr,"%g",sim->time);
1150 scmdappenddata(cmd->cmds,dataid,1,1,sim->time);
1151 for(i=1;i<nspecies*ncmpt;i++)
1152 if(i%nspecies!=0) {
1153 scmdfprintf(cmd->cmds,fptr,"%,%i",ct[i]);
1154 scmdappenddata(cmd->cmds,dataid,0,1,(double)ct[i]); }
1155 scmdfprintf(cmd->cmds,fptr,"\n");
1156 scmdflush(fptr);
1157 return CMDok;
1158
1159 scanportion:
1160 mptr=(moleculeptr) line2;
1161 for(ic=0;ic<ncmpt;ic++)
1162 if(posincompart(sim,mptr->pos,cmptlist[ic],0)) ct[ic*nspecies+mptr->ident]++;
1163 return CMDok; }
1164
1165
1166 /* cmdmolcountincmpt2 */
cmdmolcountincmpt2(simptr sim,cmdptr cmd,char * line2)1167 enum CMDcode cmdmolcountincmpt2(simptr sim,cmdptr cmd,char *line2) {
1168 FILE *fptr;
1169 char nm[STRCHAR],state[STRCHAR];
1170 compartssptr cmptss;
1171 int itct,c,i,nspecies,er,dataid;
1172 moleculeptr mptr;
1173 enum MolecState ms;
1174 static compartptr cmpt;
1175 static int inscan=0,*ct;
1176
1177 if(inscan) goto scanportion;
1178 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
1179
1180 SCMDCHECK(cmd->i1!=-1,"error on setup"); // failed before, don't try again
1181 cmptss=sim->cmptss;
1182 SCMDCHECK(cmptss,"no compartments defined");
1183 SCMDCHECK(sim->mols,"molecules are undefined");
1184 SCMDCHECK(line2,"missing argument");
1185 itct=sscanf(line2,"%s %s",nm,state);
1186 SCMDCHECK(itct==2,"cannot read arguments");
1187 c=stringfind(cmptss->cnames,cmptss->ncmpt,nm);
1188 SCMDCHECK(c>=0,"compartment name not recognized");
1189 ms=molstring2ms(state);
1190 SCMDCHECK(ms!=MSnone,"molecule state not recognized");
1191 SCMDCHECK(ms!=MSbsoln,"bsoln molecule state not permitted");
1192 cmpt=cmptss->cmptlist[c];
1193 line2=strnword(line2,3);
1194 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
1195 SCMDCHECK(er!=-1,"file or data name not recognized");
1196
1197 nspecies=sim->mols->nspecies;
1198 if(cmd->i1!=nspecies) { // allocate counter if required
1199 cmdv1free(cmd);
1200 cmd->i1=nspecies;
1201 cmd->freefn=&cmdv1free;
1202 cmd->v1=calloc(nspecies,sizeof(int));
1203 if(!cmd->v1) {cmd->i1=-1;return CMDwarn;} }
1204
1205 ct=(int*)cmd->v1;
1206 for(i=0;i<nspecies;i++) ct[i]=0;
1207
1208 inscan=1;
1209 molscancmd(sim,-1,NULL,ms,cmd,cmdmolcountincmpt2);
1210 inscan=0;
1211
1212 scmdfprintf(cmd->cmds,fptr,"%g",sim->time);
1213 scmdappenddata(cmd->cmds,dataid,1,1,sim->time);
1214 for(i=1;i<nspecies;i++) {
1215 scmdfprintf(cmd->cmds,fptr,"%,%i",ct[i]);
1216 scmdappenddata(cmd->cmds,dataid,0,1,(double)ct[i]); }
1217 scmdfprintf(cmd->cmds,fptr,"\n");
1218 scmdflush(fptr);
1219 return CMDok;
1220
1221 scanportion:
1222 mptr=(moleculeptr) line2;
1223 if(posincompart(sim,mptr->pos,cmpt,0)) ct[mptr->ident]++;
1224 return CMDok; }
1225
1226
1227 /* cmdmolcountonsurf */
cmdmolcountonsurf(simptr sim,cmdptr cmd,char * line2)1228 enum CMDcode cmdmolcountonsurf(simptr sim,cmdptr cmd,char *line2) {
1229 FILE *fptr;
1230 char nm[STRCHAR];
1231 surfacessptr srfss;
1232 int itct,s,i,nspecies,er,dataid;
1233 moleculeptr mptr;
1234 static int inscan=0,*ct;
1235 static surfaceptr srf;
1236
1237 if(inscan) goto scanportion;
1238 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
1239
1240 SCMDCHECK(cmd->i1!=-1,"error on setup"); // failed before, don't try again
1241 srfss=sim->srfss;
1242 SCMDCHECK(srfss,"no surfaces defined");
1243 SCMDCHECK(sim->mols,"molecules are undefined");
1244 SCMDCHECK(line2,"missing argument");
1245 itct=sscanf(line2,"%s",nm);
1246 SCMDCHECK(itct==1,"cannot read argument");
1247 s=stringfind(srfss->snames,srfss->nsrf,nm);
1248 SCMDCHECK(s>=0,"surface name '%s' not recognized",nm);
1249 srf=srfss->srflist[s];
1250 line2=strnword(line2,2);
1251 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
1252 SCMDCHECK(er!=-1,"file or data name not recognized");
1253
1254 nspecies=sim->mols->nspecies;
1255 if(cmd->i1!=nspecies) { // allocate counter if required
1256 cmdv1free(cmd);
1257 cmd->i1=nspecies;
1258 cmd->freefn=&cmdv1free;
1259 cmd->v1=calloc(nspecies,sizeof(int));
1260 if(!cmd->v1) {cmd->i1=-1;return CMDwarn;} }
1261
1262 ct=(int*)cmd->v1;
1263 for(i=0;i<nspecies;i++) ct[i]=0;
1264
1265 inscan=1;
1266 molscancmd(sim,-1,NULL,MSall,cmd,cmdmolcountonsurf);
1267 inscan=0;
1268
1269 scmdfprintf(cmd->cmds,fptr,"%g",sim->time);
1270 scmdappenddata(cmd->cmds,dataid,1,1,sim->time);
1271 for(i=1;i<nspecies;i++) {
1272 scmdfprintf(cmd->cmds,fptr,"%,%i",ct[i]);
1273 scmdappenddata(cmd->cmds,dataid,0,1,(double)ct[i]); }
1274 scmdfprintf(cmd->cmds,fptr,"\n");
1275 scmdflush(fptr);
1276 return CMDok;
1277
1278 scanportion:
1279 mptr=(moleculeptr) line2;
1280 if(mptr->mstate!=MSsoln && mptr->pnl->srf==srf) ct[mptr->ident]++;
1281 return CMDok; }
1282
1283
1284 /* cmdmolcountspace */
cmdmolcountspace(simptr sim,cmdptr cmd,char * line2)1285 enum CMDcode cmdmolcountspace(simptr sim,cmdptr cmd,char *line2) {
1286 FILE *fptr;
1287 int dim,i,itct,ax2,d,bin,average,*ctlat,ilat,*index,j,er,dataid;
1288 enum MolecState ms;
1289 char axisstr[STRCHAR];
1290 moleculeptr mptr;
1291 latticeptr lat;
1292 static double low[DIMMAX],high[DIMMAX],scale;
1293 static int inscan=0,nbin,*ct,axis;
1294
1295 if(inscan) goto scanportion;
1296 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
1297 SCMDCHECK(cmd->i1!=-1,"error on setup"); // failed before, don't try again
1298
1299 dim=sim->dim;
1300 SCMDCHECK(line2,"missing arguments");
1301 i=molstring2index1(sim,line2,&ms,&index);
1302 SCMDCHECK(i!=-1,"species is missing or cannot be read");
1303 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
1304 SCMDCHECK(i!=-3,"cannot read molecule state value");
1305 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
1306 SCMDCHECK(i!=-7,"error allocating memory");
1307 line2=strnword(line2,2);
1308 SCMDCHECK(line2,"missing arguments");
1309 itct=sscanf(line2,"%s",axisstr);
1310 SCMDCHECK(itct==1,"cannot read axis value");
1311 if(!strcmp(axisstr,"0") || !strcmp(axisstr,"x")) axis=0;
1312 else if(!strcmp(axisstr,"1") || !strcmp(axisstr,"y")) axis=1;
1313 else if(!strcmp(axisstr,"2") || !strcmp(axisstr,"z")) axis=2;
1314 else axis=3;
1315 SCMDCHECK(axis>=0 && axis<dim,"illegal axis value");
1316 line2=strnword(line2,2);
1317 SCMDCHECK(line2,"missing arguments");
1318 itct=strmathsscanf(line2,"%mlg %mlg %mi",Varnames,Varvalues,Nvar,&low[axis],&high[axis],&nbin);
1319 SCMDCHECK(itct==3,"cannot read arguments: low high bins");
1320 SCMDCHECK(low[axis]<high[axis],"low value needs to be less than high value");
1321 SCMDCHECK(nbin>0,"bins value needs to be > 0");
1322 line2=strnword(line2,4);
1323 ax2=0;
1324 for(d=0;d<dim-1;d++) {
1325 if(ax2==axis) ax2++;
1326 SCMDCHECK(line2,"missing position arguments");
1327 itct=strmathsscanf(line2,"%mlg %mlg",Varnames,Varvalues,Nvar,&low[ax2],&high[ax2]);
1328 SCMDCHECK(itct==2,"cannot read (or insufficient) position arguments");
1329 SCMDCHECK(low[ax2]<=high[ax2],"low value needs to be less than or equal to high value");
1330 line2=strnword(line2,3);
1331 ax2++; }
1332 SCMDCHECK(line2,"missing arguments");
1333 itct=strmathsscanf(line2,"%mi",Varnames,Varvalues,Nvar,&average);
1334 SCMDCHECK(itct==1,"cannot read average number");
1335 SCMDCHECK(average>=0,"illegal average value");
1336 line2=strnword(line2,2);
1337 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
1338 SCMDCHECK(er!=-1,"file or data name not recognized");
1339
1340 if(cmd->i1!=nbin) { // allocate counter if required
1341 cmdv1free(cmd);
1342 cmd->i1=nbin;
1343 cmd->freefn=&cmdv1v2free;
1344 cmd->v1=calloc(nbin,sizeof(int));
1345 if(!cmd->v1) {cmd->i1=-1;return CMDwarn;} }
1346
1347 ct=(int*)cmd->v1;
1348 if(average<=1 || cmd->invoke%average==1)
1349 for(bin=0;bin<nbin;bin++) ct[bin]=0;
1350 scale=(double)nbin/(high[axis]-low[axis]);
1351
1352 if(i!=-4) {
1353 inscan=1;
1354 molscancmd(sim,i,index,ms,cmd,cmdmolcountspace);
1355 inscan=0;
1356
1357 if(sim->latticess) {
1358 if(cmd->i2!=nbin) {
1359 free(cmd->v2);
1360 cmd->i2=nbin;
1361 cmd->v2=calloc(nbin,sizeof(int));
1362 if(!cmd->v2) {cmd->i2=-1;return CMDwarn;} }
1363
1364 ctlat=(int*)cmd->v2;
1365 for(ilat=0;ilat<sim->latticess->nlattice;ilat++) {
1366 lat=sim->latticess->latticelist[ilat];
1367 if(lat->type==LATTICEnsv) {
1368 for(j=0;j<index[PDnresults];j++) {
1369 NSV_CALL(nsv_molcountspace(lat->nsv,index[PDMAX+j],low,high,dim,nbin,axis,ctlat));
1370 for(bin=0;bin<nbin;++bin)
1371 ct[bin]+=ctlat[bin]; }}
1372 else if(lat->type==LATTICEpde) { //not implemented
1373 }}}}
1374
1375 if(average<=1) {
1376 scmdfprintf(cmd->cmds,fptr,"%g",sim->time);
1377 scmdappenddata(cmd->cmds,dataid,1,1,sim->time);
1378 for(bin=0;bin<nbin;bin++) {
1379 scmdfprintf(cmd->cmds,fptr,"%,%i",ct[bin]);
1380 scmdappenddata(cmd->cmds,dataid,0,1,(double)ct[bin]); }
1381 scmdfprintf(cmd->cmds,fptr,"\n"); }
1382 else if(cmd->invoke%average==0) {
1383 scmdfprintf(cmd->cmds,fptr,"%g",sim->time);
1384 scmdappenddata(cmd->cmds,dataid,1,1,sim->time);
1385 for(bin=0;bin<nbin;bin++) {
1386 scmdfprintf(cmd->cmds,fptr,"%,%g",(double)(ct[bin])/(double)average);
1387 scmdappenddata(cmd->cmds,dataid,0,1,(double)(ct[bin])/(double)average); }
1388 scmdfprintf(cmd->cmds,fptr,"\n"); }
1389 scmdflush(fptr);
1390 return CMDok;
1391
1392 scanportion:
1393 mptr=(moleculeptr) line2;
1394 for(d=0;d<sim->dim;d++)
1395 if(mptr->pos[d]<=low[d] || mptr->pos[d]>=high[d]) return CMDok;
1396 bin=(int)floor(scale*(mptr->pos[axis]-low[axis]));
1397 if(bin==nbin) bin--;
1398 ct[bin]++;
1399 return CMDok; }
1400
1401
1402 /* cmdmolcountspace2d */
cmdmolcountspace2d(simptr sim,cmdptr cmd,char * line2)1403 enum CMDcode cmdmolcountspace2d(simptr sim,cmdptr cmd,char *line2) {
1404 FILE *fptr;
1405 int dim,i,itct,d,bin,average,*index,curaxis,bin1,bin2,er,dataid;
1406 enum MolecState ms;
1407 char axisstr[STRCHAR];
1408 moleculeptr mptr;
1409
1410 static double low[DIMMAX],high[DIMMAX],scale1,scale2;
1411 static int inscan=0,nbin1,nbin2,*ct,axis,axis1,axis2;
1412
1413 if(inscan) goto scanportion;
1414 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
1415 SCMDCHECK(cmd->i1!=-1,"error on setup"); // failed before, don't try again
1416
1417 dim=sim->dim;
1418 SCMDCHECK(line2,"missing arguments");
1419 i=molstring2index1(sim,line2,&ms,&index);
1420 SCMDCHECK(i!=-1,"species is missing or cannot be read");
1421 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
1422 SCMDCHECK(i!=-3,"cannot read molecule state value");
1423 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
1424 SCMDCHECK(i!=-7,"error allocating memory");
1425 line2=strnword(line2,2);
1426
1427 SCMDCHECK(line2,"missing arguments");
1428 itct=sscanf(line2,"%s",axisstr); // perpendicular axis
1429 SCMDCHECK(itct==1,"cannot read axis value");
1430 if(!strcmp(axisstr,"0") || !strcmp(axisstr,"x")) axis=0;
1431 else if(!strcmp(axisstr,"1") || !strcmp(axisstr,"y")) axis=1;
1432 else if(!strcmp(axisstr,"2") || !strcmp(axisstr,"z")) axis=2;
1433 else axis=3;
1434 SCMDCHECK((dim==2 && axis==2) || (dim==3 && axis<3),"illegal axis value");
1435 line2=strnword(line2,2);
1436
1437 SCMDCHECK(line2,"missing arguments");
1438 curaxis=0;
1439 if(curaxis==axis) curaxis++; // first parallel axis
1440 itct=strmathsscanf(line2,"%mlg %mlg %mi",Varnames,Varvalues,Nvar,&low[curaxis],&high[curaxis],&nbin1);
1441 SCMDCHECK(itct==3,"cannot read arguments: low high bins");
1442 SCMDCHECK(low[curaxis]<high[curaxis],"low value needs to be less than high value");
1443 SCMDCHECK(nbin1>0,"bins value needs to be > 0");
1444 axis1=curaxis;
1445 line2=strnword(line2,4);
1446
1447 SCMDCHECK(line2,"missing arguments");
1448 curaxis++;
1449 if(curaxis==axis) curaxis++; // second parallel axis
1450 itct=strmathsscanf(line2,"%mlg %mlg %mi",Varnames,Varvalues,Nvar,&low[curaxis],&high[curaxis],&nbin2);
1451 SCMDCHECK(itct==3,"cannot read arguments: low high bins");
1452 SCMDCHECK(low[curaxis]<high[curaxis],"low value needs to be less than high value");
1453 SCMDCHECK(nbin2>0,"bins value needs to be > 0");
1454 axis2=curaxis;
1455 line2=strnword(line2,4);
1456
1457 if(dim==3) {
1458 curaxis=axis; // perpendicular axis
1459 itct=strmathsscanf(line2,"%mlg %mlg",Varnames,Varvalues,Nvar,&low[curaxis],&high[curaxis]);
1460 SCMDCHECK(itct==2,"cannot read (or insufficient) position arguments");
1461 SCMDCHECK(low[curaxis]<=high[curaxis],"low value needs to be less than or equal to high value");
1462 line2=strnword(line2,3); }
1463
1464 SCMDCHECK(line2,"missing arguments"); // average
1465 itct=strmathsscanf(line2,"%mi",Varnames,Varvalues,Nvar,&average);
1466 SCMDCHECK(itct==1,"cannot read average number");
1467 SCMDCHECK(average>=0,"illegal average value");
1468 line2=strnword(line2,2);
1469 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
1470 SCMDCHECK(er!=-1,"file or data name not recognized");
1471
1472 if(cmd->i1!=nbin1*nbin2) { // allocate counter if required
1473 cmdv1free(cmd);
1474 cmd->i1=nbin1*nbin2;
1475 cmd->freefn=&cmdv1v2free;
1476 cmd->v1=calloc(nbin1*nbin2,sizeof(int));
1477 if(!cmd->v1) {cmd->i1=-1;return CMDwarn;} }
1478
1479 ct=(int*)cmd->v1;
1480 if(average<=1 || cmd->invoke%average==1)
1481 for(bin=0;bin<nbin1*nbin2;bin++) ct[bin]=0;
1482 scale1=(double)nbin1/(high[axis1]-low[axis1]);
1483 scale2=(double)nbin2/(high[axis2]-low[axis2]);
1484
1485 if(i!=-4) {
1486 inscan=1;
1487 molscancmd(sim,i,index,ms,cmd,cmdmolcountspace2d);
1488 inscan=0; }
1489
1490 if(average<=1 || cmd->invoke%average==0) {
1491 scmdfprintf(cmd->cmds,fptr,"%g\n",sim->time);
1492 scmdappenddata(cmd->cmds,dataid,1,1,sim->time);
1493 for(bin2=0;bin2<nbin2;bin2++) {
1494 bin1=0;
1495 if(average<=1) {
1496 scmdfprintf(cmd->cmds,fptr,"%i",ct[bin2*nbin1+bin1]);
1497 scmdappenddata(cmd->cmds,dataid,0,1,(double)ct[bin2*nbin1+bin1]); }
1498 else {
1499 scmdfprintf(cmd->cmds,fptr,"%g",(double)(ct[bin2*nbin1+bin1]/(double)average));
1500 scmdappenddata(cmd->cmds,dataid,0,1,(double)(ct[bin2*nbin1+bin1]/(double)average)); }
1501 for(bin1=1;bin1<nbin1;bin1++) {
1502 if(average<=1) {
1503 scmdfprintf(cmd->cmds,fptr,"%,%i",ct[bin2*nbin1+bin1]);
1504 scmdappenddata(cmd->cmds,dataid,0,1,(double)ct[bin2*nbin1+bin1]); }
1505 else {
1506 scmdfprintf(cmd->cmds,fptr,"%,%g",(double)(ct[bin2*nbin1+bin1]/(double)average));
1507 scmdappenddata(cmd->cmds,dataid,0,1,(double)(ct[bin2*nbin1+bin1]/(double)average)); }}
1508 scmdfprintf(cmd->cmds,fptr,"\n"); }
1509 scmdflush(fptr); }
1510
1511 return CMDok;
1512
1513 scanportion:
1514 mptr=(moleculeptr) line2;
1515 for(d=0;d<sim->dim;d++)
1516 if(mptr->pos[d]<=low[d] || mptr->pos[d]>=high[d]) return CMDok;
1517 bin1=(int)floor(scale1*(mptr->pos[axis1]-low[axis1]));
1518 bin2=(int)floor(scale2*(mptr->pos[axis2]-low[axis2]));
1519 if(bin1==nbin1) bin1--;
1520 if(bin2==nbin2) bin2--;
1521 bin=bin2*nbin1+bin1;
1522 ct[bin]++;
1523 return CMDok; }
1524
1525
1526 /* cmdmolcountspaceradial */
cmdmolcountspaceradial(simptr sim,cmdptr cmd,char * line2)1527 enum CMDcode cmdmolcountspaceradial(simptr sim,cmdptr cmd,char *line2) {
1528 FILE *fptr;
1529 int i,itct,d,bin,average,*index,er,dataid;
1530 enum MolecState ms;
1531 double radius,molrad;
1532 moleculeptr mptr;
1533 static double center[DIMMAX],scale,radius2;
1534 static int inscan=0,nbin,*ct;
1535
1536 if(inscan) goto scanportion;
1537 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
1538 SCMDCHECK(cmd->i1!=-1,"error on setup"); // failed before, don't try again
1539
1540 SCMDCHECK(line2,"missing arguments");
1541 i=molstring2index1(sim,line2,&ms,&index);
1542 SCMDCHECK(i!=-1,"species is missing or cannot be read");
1543 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
1544 SCMDCHECK(i!=-3,"cannot read molecule state value");
1545 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
1546 SCMDCHECK(i!=-7,"error allocating memory");
1547 line2=strnword(line2,2);
1548 SCMDCHECK(line2,"missing arguments");
1549 for(d=0;d<sim->dim;d++) {
1550 itct=strmathsscanf(line2,"%mlg",Varnames,Varvalues,Nvar,¢er[d]);
1551 SCMDCHECK(itct==1,"missing center value");
1552 line2=strnword(line2,2);
1553 SCMDCHECK(line2,"missing arguments"); }
1554 itct=strmathsscanf(line2,"%mlg %mi",Varnames,Varvalues,Nvar,&radius,&nbin);
1555 SCMDCHECK(itct==2,"cannot read arguments: radius bins");
1556 SCMDCHECK(radius>0,"radius needs to be greater than 0");
1557 SCMDCHECK(nbin>0,"bins value needs to be > 0");
1558 line2=strnword(line2,3);
1559 SCMDCHECK(line2,"missing arguments");
1560 itct=strmathsscanf(line2,"%mi",Varnames,Varvalues,Nvar,&average);
1561 SCMDCHECK(itct==1,"cannot read average number");
1562 SCMDCHECK(average>=0,"illegal average value");
1563 line2=strnword(line2,2);
1564 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
1565 SCMDCHECK(er!=-1,"file or data name not recognized");
1566
1567 if(cmd->i1!=nbin) { // allocate counter if required
1568 cmdv1free(cmd);
1569 cmd->i1=nbin;
1570 cmd->freefn=&cmdv1v2free;
1571 cmd->v1=calloc(nbin,sizeof(int));
1572 if(!cmd->v1) {cmd->i1=-1;return CMDwarn;} }
1573
1574 ct=(int*)cmd->v1;
1575 if(average<=1 || cmd->invoke%average==1)
1576 for(bin=0;bin<nbin;bin++) ct[bin]=0;
1577 scale=(double)nbin/radius;
1578
1579 radius2=radius*radius;
1580
1581 if(i!=-4) {
1582 inscan=1;
1583 molscancmd(sim,i,index,ms,cmd,cmdmolcountspaceradial);
1584 inscan=0; }
1585
1586 if(average<=1) {
1587 scmdfprintf(cmd->cmds,fptr,"%g",sim->time);
1588 scmdappenddata(cmd->cmds,dataid,1,1,sim->time);
1589 for(bin=0;bin<nbin;bin++) {
1590 scmdfprintf(cmd->cmds,fptr,"%,%i",ct[bin]);
1591 scmdappenddata(cmd->cmds,dataid,0,1,(double)ct[bin]); }
1592 scmdfprintf(cmd->cmds,fptr,"\n"); }
1593 else if(cmd->invoke%average==0) {
1594 scmdfprintf(cmd->cmds,fptr,"%g",sim->time);
1595 scmdappenddata(cmd->cmds,dataid,1,1,sim->time);
1596 for(bin=0;bin<nbin;bin++) {
1597 scmdfprintf(cmd->cmds,fptr,"%,%g",(double)(ct[bin])/(double)average);
1598 scmdappenddata(cmd->cmds,dataid,0,1,(double)(ct[bin])/(double)average); }
1599 scmdfprintf(cmd->cmds,fptr,"\n"); }
1600 scmdflush(fptr);
1601 return CMDok;
1602
1603 scanportion:
1604 mptr=(moleculeptr) line2;
1605 molrad=0;
1606 for(d=0;d<sim->dim;d++)
1607 molrad+=(mptr->pos[d]-center[d])*(mptr->pos[d]-center[d]);
1608 if(molrad<radius2) {
1609 molrad=sqrt(molrad);
1610 bin=(int)floor(scale*molrad);
1611 if(bin==nbin) bin--;
1612 ct[bin]++; }
1613 return CMDok; }
1614
1615
1616 /* cmdmolcountspacepolarangle */
cmdmolcountspacepolarangle(simptr sim,cmdptr cmd,char * line2)1617 enum CMDcode cmdmolcountspacepolarangle(simptr sim,cmdptr cmd,char *line2) {
1618 FILE *fptr;
1619 int i,itct,d,bin,average,*index,er,dataid;
1620 enum MolecState ms;
1621 double radiusmin,radiusmax,molrad,poleleninv,angle;
1622 moleculeptr mptr;
1623 static double center[DIMMAX],pole[DIMMAX],poleangle,scale,radiusmin2,radiusmax2;
1624 static int inscan=0,nbin,*ct;
1625
1626 if(inscan) goto scanportion;
1627 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
1628 SCMDCHECK(cmd->i1!=-1,"error on setup"); // failed before, don't try again
1629
1630 SCMDCHECK(line2,"missing arguments");
1631 i=molstring2index1(sim,line2,&ms,&index);
1632 SCMDCHECK(i!=-1,"species is missing or cannot be read");
1633 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
1634 SCMDCHECK(i!=-3,"cannot read molecule state value");
1635 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
1636 SCMDCHECK(i!=-7,"error allocating memory");
1637 line2=strnword(line2,2);
1638 SCMDCHECK(line2,"missing arguments");
1639 for(d=0;d<sim->dim;d++) {
1640 itct=strmathsscanf(line2,"%mlg",Varnames,Varvalues,Nvar,¢er[d]);
1641 SCMDCHECK(itct==1,"missing center value");
1642 line2=strnword(line2,2);
1643 SCMDCHECK(line2,"missing arguments"); }
1644 for(d=0;d<sim->dim;d++) {
1645 itct=strmathsscanf(line2,"%mlg",Varnames,Varvalues,Nvar,&pole[d]);
1646 SCMDCHECK(itct==1,"missing pole value");
1647 line2=strnword(line2,2);
1648 SCMDCHECK(line2,"missing arguments"); }
1649 itct=strmathsscanf(line2,"%mlg %mlg %mi",Varnames,Varvalues,Nvar,&radiusmin,&radiusmax,&nbin);
1650 SCMDCHECK(itct==3,"cannot read arguments: radius_min radius_max bins");
1651 SCMDCHECK(nbin>0,"bins value needs to be > 0");
1652 line2=strnword(line2,4);
1653 SCMDCHECK(line2,"missing arguments");
1654 itct=strmathsscanf(line2,"%mi",Varnames,Varvalues,Nvar,&average);
1655 SCMDCHECK(itct==1,"cannot read average number");
1656 SCMDCHECK(average>=0,"illegal average value");
1657 line2=strnword(line2,2);
1658 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
1659 SCMDCHECK(er!=-1,"file or data name not recognized");
1660
1661 if(cmd->i1!=nbin) { // allocate counter if required
1662 cmdv1free(cmd);
1663 cmd->i1=nbin;
1664 cmd->freefn=&cmdv1v2free;
1665 cmd->v1=calloc(nbin,sizeof(int));
1666 if(!cmd->v1) {cmd->i1=-1;return CMDwarn;} }
1667
1668 ct=(int*)cmd->v1;
1669 if(average<=1 || cmd->invoke%average==1)
1670 for(bin=0;bin<nbin;bin++) ct[bin]=0;
1671 scale=(double)nbin/(sim->dim==2?2*PI:PI);
1672
1673 radiusmin2=radiusmin>=0?radiusmin*radiusmin:0;
1674 radiusmax2=radiusmax>=0?radiusmax*radiusmax:-1;
1675 if(sim->dim==2) {
1676 SCMDCHECK(pole[0]!=0 || pole[1]!=0,"pole vector is equal to zero");
1677 poleangle=atan2(pole[1],pole[0]); }
1678 else {
1679 poleangle=0;
1680 poleleninv=sqrt(pole[0]*pole[0]+pole[1]*pole[1]+pole[2]*pole[2]);
1681 SCMDCHECK(poleleninv>0,"pole vector is equal to zero");
1682 poleleninv=1.0/poleleninv;
1683 pole[0]*=poleleninv;
1684 pole[1]*=poleleninv;
1685 pole[2]*=poleleninv; }
1686
1687 if(i!=-4) {
1688 inscan=1;
1689 molscancmd(sim,i,index,ms,cmd,cmdmolcountspacepolarangle);
1690 inscan=0; }
1691
1692 if(average<=1) {
1693 scmdfprintf(cmd->cmds,fptr,"%g",sim->time);
1694 scmdappenddata(cmd->cmds,dataid,1,1,sim->time);
1695 for(bin=0;bin<nbin;bin++) {
1696 scmdfprintf(cmd->cmds,fptr,"%,%i",ct[bin]);
1697 scmdappenddata(cmd->cmds,dataid,0,1,(double)ct[bin]); }
1698 scmdfprintf(cmd->cmds,fptr,"\n"); }
1699 else if(cmd->invoke%average==0) {
1700 scmdfprintf(cmd->cmds,fptr,"%g",sim->time);
1701 scmdappenddata(cmd->cmds,dataid,1,1,sim->time);
1702 for(bin=0;bin<nbin;bin++) {
1703 scmdfprintf(cmd->cmds,fptr,"%,%g",(double)(ct[bin])/(double)average);
1704 scmdappenddata(cmd->cmds,dataid,0,1,(double)(ct[bin])/(double)average); }
1705 scmdfprintf(cmd->cmds,fptr,"\n"); }
1706 scmdflush(fptr);
1707 return CMDok;
1708
1709 scanportion:
1710 mptr=(moleculeptr) line2;
1711 molrad=0;
1712 for(d=0;d<sim->dim;d++)
1713 molrad+=(mptr->pos[d]-center[d])*(mptr->pos[d]-center[d]);
1714 if(molrad>=radiusmin2 && (radiusmax2==-1 || molrad<=radiusmax2)) {
1715 if(sim->dim==2) {
1716 angle=atan2(mptr->pos[1]-center[1],mptr->pos[0]-center[0]);
1717 angle-=poleangle;
1718 if(angle<0) angle+=2*PI;
1719 else if(angle>2*PI) angle-=2*PI; }
1720 else {
1721 angle=acos(((mptr->pos[0]-center[0])*pole[0]+(mptr->pos[1]-center[1])*pole[1]+(mptr->pos[2]-center[2])*pole[2])/sqrt(molrad)); }
1722 bin=(int)floor(scale*angle);
1723 if(bin==nbin) bin--;
1724 ct[bin]++; }
1725 return CMDok; }
1726
1727
1728 /* cmdradialdistribution */
cmdradialdistribution(simptr sim,cmdptr cmd,char * line2)1729 enum CMDcode cmdradialdistribution(simptr sim,cmdptr cmd,char *line2) {
1730 FILE *fptr;
1731 int dim,i1,itct,d,bin,average,*index1,i,ll,m,wrap[DIMMAX],er,dataid;
1732 enum MolecState ms1,mslo,mshi,ms;
1733 moleculeptr mptr,mptr2;
1734 boxptr bptr;
1735 double dist,scale2,rdf;
1736 static double scale,radius,syswidth[DIMMAX];
1737 static int inscan=0,nbin,*ct,i2,*index2,lllo,llhi,mcount;
1738 static enum MolecState ms2;
1739
1740 if(inscan) goto scanportion;
1741 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
1742 SCMDCHECK(cmd->i1!=-1,"error on setup"); // failed before, don't try again
1743
1744 SCMDCHECK(line2,"missing arguments");
1745 i1=molstring2index1(sim,line2,&ms1,&index1);
1746 SCMDCHECK(i1!=-1,"species is missing or cannot be read");
1747 SCMDCHECK(i1!=-2,"mismatched or improper parentheses around molecule state");
1748 SCMDCHECK(i1!=-3,"cannot read molecule state value");
1749 SCMDCHECK(i1!=-4,"molecule name not recognized");
1750 SCMDCHECK(i1!=-7,"error allocating memory");
1751 line2=strnword(line2,2);
1752 SCMDCHECK(line2,"missing arguments");
1753 i2=molstring2index1(sim,line2,&ms2,&index2);
1754 SCMDCHECK(i2!=-1,"species is missing or cannot be read");
1755 SCMDCHECK(i2!=-2,"mismatched or improper parentheses around molecule state");
1756 SCMDCHECK(i2!=-3,"cannot read molecule state value");
1757 SCMDCHECK(i2!=-4,"molecule name not recognized");
1758 SCMDCHECK(i2!=-7,"error allocating memory");
1759 line2=strnword(line2,2);
1760 SCMDCHECK(line2,"missing arguments");
1761 itct=strmathsscanf(line2,"%mlg %mi %mi",Varnames,Varvalues,Nvar,&radius,&nbin,&average);
1762 SCMDCHECK(itct==3,"cannot read arguments: radius bins average");
1763 SCMDCHECK(radius>0,"radius needs to be greater than 0");
1764 SCMDCHECK(nbin>0,"bins value needs to be > 0");
1765 SCMDCHECK(average>=0,"illegal average value");
1766 line2=strnword(line2,4);
1767 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
1768 SCMDCHECK(er!=-1,"file or data name not recognized");
1769
1770 if(cmd->i1!=nbin) { // allocate counter if required
1771 cmdv1free(cmd);
1772 cmd->i1=nbin;
1773 cmd->freefn=&cmdv1v2free;
1774 cmd->v1=calloc(nbin,sizeof(int));
1775 if(!cmd->v1) {cmd->i1=-1;return CMDwarn;} }
1776
1777 dim=sim->dim;
1778 ct=(int*)cmd->v1;
1779 if(average<=1 || cmd->invoke%average==1) {
1780 for(bin=0;bin<nbin;bin++) ct[bin]=0;
1781 mcount=0; }
1782 scale=(double)nbin/radius; // 1/scale is radial distance per bin
1783 if(dim==1) scale2=2.0/scale; // scale2*bin^dim = volume inside bin where bin=0 for first bin
1784 else if(dim==2) scale2=PI/(scale*scale);
1785 else scale2=(4.0*PI/3.0)/(scale*scale*scale);
1786
1787 lllo=llhi=-1;
1788 if(ms2<MSMAX) {
1789 mslo=ms2;
1790 mshi=(enum MolecState)(ms2+1); }
1791 else {
1792 mslo=(enum MolecState) 0;
1793 mshi=(enum MolecState) MSMAX; }
1794 for(i=0;i<index2[PDnresults];i++)
1795 for(ms=mslo;ms<mshi;ms=(enum MolecState)(ms+1)) {
1796 ll=sim->mols->listlookup[index2[PDMAX+i]][ms];
1797 if(ll<lllo || lllo==-1) lllo=ll;
1798 if(ll>=llhi || llhi==-1) llhi=ll+1; }
1799
1800 for(d=0;d<dim;d++)
1801 syswidth[d]=sim->wlist[2*d+1]->pos-sim->wlist[2*d]->pos;
1802
1803 inscan=1;
1804 molscancmd(sim,i1,index1,ms1,cmd,cmdradialdistribution);
1805 inscan=0;
1806
1807 if(average<=1 || cmd->invoke%average==0) {
1808 if(average<1) average=1;
1809 scmdfprintf(cmd->cmds,fptr,"%g",sim->time);
1810 scmdappenddata(cmd->cmds,dataid,1,1,sim->time);
1811 for(bin=0;bin<nbin;bin++) {
1812 if(dim==1) rdf=(double)ct[bin]/((double)mcount*scale2);
1813 else if(dim==2) rdf=(double)ct[bin]/((double)mcount*scale2*(2*bin+1)); // (bin+1)^2-bin^2=(2*bin+1)
1814 else rdf=(double)ct[bin]/((double)mcount*scale2*(3*bin*bin+3*bin+1)); // (bin+1)^3-bin^3=(3*bin^2+3*bin+1)
1815 scmdfprintf(cmd->cmds,fptr,"%,%g",rdf);
1816 scmdappenddata(cmd->cmds,dataid,0,1,rdf); }
1817 scmdfprintf(cmd->cmds,fptr,"\n"); }
1818 scmdflush(fptr);
1819 return CMDok;
1820
1821 scanportion:
1822 dim=sim->dim;
1823 mptr=(moleculeptr) line2;
1824 mcount++; // mcount is number of "center" molecules
1825 bptr=boxscansphere(sim,mptr->pos,radius,NULL,wrap);
1826 while(bptr) {
1827 for(ll=lllo;ll<llhi;ll++)
1828 for(m=0;m<bptr->nmol[ll];m++) {
1829 mptr2=bptr->mol[ll][m];
1830 if(mptr2!=mptr && molismatch(mptr2,i2,index2,ms2)) {
1831 dist=0;
1832 for(d=0;d<dim;d++)
1833 dist+=(mptr2->pos[d]+wrap[d]*syswidth[d]-mptr->pos[d])*(mptr2->pos[d]+wrap[d]*syswidth[d]-mptr->pos[d]);
1834 dist=sqrt(dist);
1835 bin=(int)floor(scale*dist);
1836 if(bin<nbin)
1837 ct[bin]++; }}
1838 bptr=boxscansphere(sim,mptr->pos,radius,bptr,wrap); }
1839 return CMDok; }
1840
1841
1842 /* cmdradialdistribution2 */
cmdradialdistribution2(simptr sim,cmdptr cmd,char * line2)1843 enum CMDcode cmdradialdistribution2(simptr sim,cmdptr cmd,char *line2) {
1844 FILE *fptr;
1845 int dim,i1,itct,d,bin,average,*index1,i,ll,m,wrap[DIMMAX],er,dataid;
1846 enum MolecState ms1,mslo,mshi,ms;
1847 moleculeptr mptr,mptr2;
1848 boxptr bptr;
1849 double dist,scale2,rdf;
1850 static double scale,radius,syswidth[DIMMAX],lowpos[DIMMAX],highpos[DIMMAX];
1851 static int inscan=0,nbin,*ct,i2,*index2,lllo,llhi,mcount;
1852 static enum MolecState ms2;
1853
1854 if(inscan) goto scanportion;
1855 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
1856 SCMDCHECK(cmd->i1!=-1,"error on setup"); // failed before, don't try again
1857
1858 SCMDCHECK(line2,"missing arguments");
1859 i1=molstring2index1(sim,line2,&ms1,&index1);
1860 SCMDCHECK(i1!=-1,"species is missing or cannot be read");
1861 SCMDCHECK(i1!=-2,"mismatched or improper parentheses around molecule state");
1862 SCMDCHECK(i1!=-3,"cannot read molecule state value");
1863 SCMDCHECK(i1!=-4,"molecule name not recognized");
1864 SCMDCHECK(i1!=-7,"error allocating memory");
1865 line2=strnword(line2,2);
1866 SCMDCHECK(line2,"missing arguments");
1867 i2=molstring2index1(sim,line2,&ms2,&index2);
1868 SCMDCHECK(i2!=-1,"species is missing or cannot be read");
1869 SCMDCHECK(i2!=-2,"mismatched or improper parentheses around molecule state");
1870 SCMDCHECK(i2!=-3,"cannot read molecule state value");
1871 SCMDCHECK(i2!=-4,"molecule name not recognized");
1872 SCMDCHECK(i2!=-7,"error allocating memory");
1873 line2=strnword(line2,2);
1874 SCMDCHECK(line2,"missing arguments");
1875 for(d=0;d<sim->dim;d++) {
1876 itct=strmathsscanf(line2,"%mlg %mlg",Varnames,Varvalues,Nvar,&lowpos[d],&highpos[d]);
1877 SCMDCHECK(itct==2,"missing arguments");
1878 SCMDCHECK(lowpos[d]<=highpos[d],"low position value needs to be <= high position value");
1879 line2=strnword(line2,3);
1880 SCMDCHECK(line2,"missing arguments"); }
1881 itct=strmathsscanf(line2,"%mlg %mi %mi",Varnames,Varvalues,Nvar,&radius,&nbin,&average);
1882 SCMDCHECK(itct==3,"cannot read arguments: radius bins average");
1883 SCMDCHECK(radius>0,"radius needs to be greater than 0");
1884 SCMDCHECK(nbin>0,"bins value needs to be > 0");
1885 SCMDCHECK(average>=0,"illegal average value");
1886 line2=strnword(line2,4);
1887 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
1888 SCMDCHECK(er!=-1,"file or data name not recognized");
1889
1890 if(cmd->i1!=nbin) { // allocate counter if required
1891 cmdv1free(cmd);
1892 cmd->i1=nbin;
1893 cmd->freefn=&cmdv1v2free;
1894 cmd->v1=calloc(nbin,sizeof(int));
1895 if(!cmd->v1) {cmd->i1=-1;return CMDwarn;} }
1896
1897 dim=sim->dim;
1898 ct=(int*)cmd->v1;
1899 if(average<=1 || cmd->invoke%average==1) {
1900 for(bin=0;bin<nbin;bin++) ct[bin]=0;
1901 mcount=0; }
1902 scale=(double)nbin/radius; // 1/scale is radial distance per bin
1903 if(dim==1) scale2=2.0/scale; // scale2*bin^dim = volume inside bin where bin=0 for first bin
1904 else if(dim==2) scale2=PI/(scale*scale);
1905 else scale2=(4.0*PI/3.0)/(scale*scale*scale);
1906
1907 lllo=llhi=-1;
1908 if(ms2<MSMAX) {
1909 mslo=ms2;
1910 mshi=(enum MolecState)(ms2+1); }
1911 else {
1912 mslo=(enum MolecState) 0;
1913 mshi=(enum MolecState) MSMAX; }
1914 for(i=0;i<index2[PDnresults];i++)
1915 for(ms=mslo;ms<mshi;ms=(enum MolecState)(ms+1)) {
1916 ll=sim->mols->listlookup[index2[PDMAX+i]][ms];
1917 if(ll<lllo || lllo==-1) lllo=ll;
1918 if(ll>=llhi || llhi==-1) llhi=ll+1; }
1919
1920 for(d=0;d<dim;d++)
1921 syswidth[d]=sim->wlist[2*d+1]->pos-sim->wlist[2*d]->pos;
1922
1923 inscan=1;
1924 molscancmd(sim,i1,index1,ms1,cmd,cmdradialdistribution2);
1925 inscan=0;
1926
1927 if(average<=1 || cmd->invoke%average==0) {
1928 if(average<1) average=1;
1929 scmdfprintf(cmd->cmds,fptr,"%g",sim->time);
1930 scmdappenddata(cmd->cmds,dataid,1,1,sim->time);
1931 for(bin=0;bin<nbin;bin++) {
1932 if(dim==1) rdf=(double)ct[bin]/((double)mcount*scale2);
1933 else if(dim==2) rdf=(double)ct[bin]/((double)mcount*scale2*(2*bin+1)); // (bin+1)^2-bin^2=(2*bin+1)
1934 else rdf=(double)ct[bin]/((double)mcount*scale2*(3*bin*bin+3*bin+1)); // (bin+1)^3-bin^3=(3*bin^2+3*bin+1)
1935 scmdfprintf(cmd->cmds,fptr,"%,%g",rdf);
1936 scmdappenddata(cmd->cmds,dataid,0,1,rdf); }
1937 scmdfprintf(cmd->cmds,fptr,"\n"); }
1938 scmdflush(fptr);
1939 return CMDok;
1940
1941 scanportion:
1942 dim=sim->dim;
1943 mptr=(moleculeptr) line2;
1944 for(d=0;d<dim;d++) {
1945 if(mptr->pos[d]<lowpos[d] || mptr->pos[d]>highpos[d]) return CMDok; }
1946 mcount++; // mcount is number of "center" molecules
1947 bptr=boxscansphere(sim,mptr->pos,radius,NULL,wrap);
1948 while(bptr) {
1949 for(ll=lllo;ll<llhi;ll++)
1950 for(m=0;m<bptr->nmol[ll];m++) {
1951 mptr2=bptr->mol[ll][m];
1952 if(mptr2!=mptr && molismatch(mptr2,i2,index2,ms2)) {
1953 dist=0;
1954 for(d=0;d<dim;d++)
1955 dist+=(mptr2->pos[d]+wrap[d]*syswidth[d]-mptr->pos[d])*(mptr2->pos[d]+wrap[d]*syswidth[d]-mptr->pos[d]);
1956 dist=sqrt(dist);
1957 bin=(int)floor(scale*dist);
1958 if(bin<nbin)
1959 ct[bin]++; }}
1960 bptr=boxscansphere(sim,mptr->pos,radius,bptr,wrap); }
1961 return CMDok; }
1962
1963
1964 /* cmdmolcountspecies */
cmdmolcountspecies(simptr sim,cmdptr cmd,char * line2)1965 enum CMDcode cmdmolcountspecies(simptr sim,cmdptr cmd,char *line2) {
1966 FILE *fptr;
1967 int i,*index,count,er,dataid;
1968 enum MolecState ms;
1969
1970 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
1971 i=molstring2index1(sim,line2,&ms,&index);
1972 SCMDCHECK(i!=-1,"species is missing or cannot be read");
1973 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
1974 SCMDCHECK(i!=-3,"cannot read molecule state value");
1975 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
1976 SCMDCHECK(i!=-7,"error allocating memory");
1977 line2=strnword(line2,2);
1978 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
1979 SCMDCHECK(er!=-1,"file or data name not recognized");
1980
1981 count=(i==-4)?0:molcount(sim,i,index,ms,-1);
1982 scmdfprintf(cmd->cmds,fptr,"%g%,%i\n",sim->time,count);
1983 scmdappenddata(cmd->cmds,dataid,1,2,sim->time,(double)count);
1984 scmdflush(fptr);
1985 return CMDok; }
1986
1987
1988 /* cmdmolcountspecieslist */
cmdmolcountspecieslist(simptr sim,cmdptr cmd,char * line2)1989 enum CMDcode cmdmolcountspecieslist(simptr sim,cmdptr cmd,char *line2) {
1990 FILE *fptr;
1991 int i,*index,count,er,dataid;
1992 enum MolecState ms;
1993
1994 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
1995 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
1996 SCMDCHECK(er!=-1,"file or data name not recognized");
1997 scmdfprintf(cmd->cmds,fptr,"%g",sim->time);
1998 scmdappenddata(cmd->cmds,dataid,1,1,sim->time);
1999 while((line2=strnword(line2,2))) {
2000 i=molstring2index1(sim,line2,&ms,&index);
2001 SCMDCHECK(i!=-1,"species is missing or cannot be read");
2002 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
2003 SCMDCHECK(i!=-3,"cannot read molecule state value");
2004 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
2005 SCMDCHECK(i!=-7,"error allocating memory");
2006 count=(i==-4)?0:molcount(sim,i,index,ms,-1);
2007 scmdfprintf(cmd->cmds,fptr,"%,%i",count);
2008 scmdappenddata(cmd->cmds,dataid,0,1,(double)count); }
2009
2010 scmdfprintf(cmd->cmds,fptr,"\n");
2011 scmdflush(fptr);
2012 return CMDok; }
2013
2014
2015 /* cmdmollistsize */
cmdmollistsize(simptr sim,cmdptr cmd,char * line2)2016 enum CMDcode cmdmollistsize(simptr sim,cmdptr cmd,char *line2) {
2017 FILE *fptr;
2018 int ll,itct,er,dataid;
2019 char listname[STRCHAR];
2020
2021 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
2022 itct=sscanf(line2,"%s",listname);
2023 SCMDCHECK(itct==1,"cannot read molecule list name");
2024 SCMDCHECK(sim->mols && sim->mols->nlist>0,"no molecule lists defined");
2025 ll=stringfind(sim->mols->listname,sim->mols->nlist,listname);
2026 SCMDCHECK(ll>=0,"molecule list name not recognized");
2027 line2=strnword(line2,2);
2028 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
2029 SCMDCHECK(er!=-1,"file or data name not recognized");
2030 scmdfprintf(cmd->cmds,fptr,"%g%,%i\n",sim->time,sim->mols->nl[ll]);
2031 scmdappenddata(cmd->cmds,dataid,1,2,sim->time,(double)(sim->mols->nl[ll]));
2032 scmdflush(fptr);
2033 return CMDok; }
2034
2035
2036 /* cmdlistmols */
cmdlistmols(simptr sim,cmdptr cmd,char * line2)2037 enum CMDcode cmdlistmols(simptr sim,cmdptr cmd,char *line2) {
2038 int d,er;
2039 char string[STRCHAR];
2040 moleculeptr mptr;
2041 static FILE *fptr;
2042 static int inscan=0;
2043 static int dataid=-1;
2044
2045 if(inscan) goto scanportion;
2046 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
2047 SCMDCHECK(sim->mols,"molecules are undefined");
2048
2049 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
2050 SCMDCHECK(er!=-1,"file or data name not recognized");
2051
2052 inscan=1;
2053 molscancmd(sim,-1,NULL,MSall,cmd,cmdlistmols);
2054 inscan=0;
2055
2056 scmdflush(fptr);
2057 return CMDok;
2058
2059 scanportion:
2060 mptr=(moleculeptr) line2;
2061 scmdfprintf(cmd->cmds,fptr,"%s(%s)",sim->mols->spname[mptr->ident],molms2string(mptr->mstate,string));
2062 scmdappenddata(cmd->cmds,dataid,1,2,(double)(mptr->ident),(double)(mptr->mstate));
2063 for(d=0;d<sim->dim;d++) {
2064 scmdfprintf(cmd->cmds,fptr,"%,%g",mptr->pos[d]);
2065 scmdappenddata(cmd->cmds,dataid,0,1,mptr->pos[d]); }
2066 scmdfprintf(cmd->cmds,fptr,"%,%s\n",molserno2string(mptr->serno,string));
2067 scmdappenddata(cmd->cmds,dataid,0,1,(double)(mptr->serno));
2068 return CMDok; }
2069
2070
2071 /* cmdlistmols2 */
cmdlistmols2(simptr sim,cmdptr cmd,char * line2)2072 enum CMDcode cmdlistmols2(simptr sim,cmdptr cmd,char *line2) {
2073 int d,er;
2074 moleculeptr mptr;
2075 static FILE *fptr;
2076 static int inscan=0,invk,dataid=-1;
2077 char string[STRCHAR];
2078
2079 if(inscan) goto scanportion;
2080 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
2081 SCMDCHECK(sim->mols,"molecules are undefined");
2082
2083 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
2084 SCMDCHECK(er!=-1,"file or data name not recognized");
2085 invk=cmd?cmd->invoke:0;
2086
2087 inscan=1;
2088 molscancmd(sim,-1,NULL,MSall,cmd,cmdlistmols2);
2089 inscan=0;
2090
2091 scmdflush(fptr);
2092 return CMDok;
2093
2094 scanportion:
2095 mptr=(moleculeptr) line2;
2096 scmdfprintf(cmd->cmds,fptr,"%i%,%i%,%i",invk,mptr->ident,mptr->mstate);
2097 scmdappenddata(cmd->cmds,dataid,1,3,(double)invk,(double)(mptr->ident),(double)(mptr->mstate));
2098 for(d=0;d<sim->dim;d++) {
2099 scmdfprintf(cmd->cmds,fptr,"%,%g",mptr->pos[d]);
2100 scmdappenddata(cmd->cmds,dataid,0,1,mptr->pos[d]); }
2101 scmdfprintf(cmd->cmds,fptr,"%,%s\n",molserno2string(mptr->serno,string));
2102 scmdappenddata(cmd->cmds,dataid,0,1,(double)(mptr->serno));
2103 return CMDok; }
2104
2105
2106 /* cmdlistmols3 */
cmdlistmols3(simptr sim,cmdptr cmd,char * line2)2107 enum CMDcode cmdlistmols3(simptr sim,cmdptr cmd,char *line2) {
2108 int i,*index,d,er;
2109 moleculeptr mptr;
2110 enum MolecState ms;
2111 static FILE *fptr;
2112 static int inscan=0,invk,dataid=-1;
2113 char string[STRCHAR];
2114
2115 if(inscan) goto scanportion;
2116 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
2117
2118 i=molstring2index1(sim,line2,&ms,&index);
2119 SCMDCHECK(i!=-1,"species is missing or cannot be read");
2120 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
2121 SCMDCHECK(i!=-3,"cannot read molecule state value");
2122 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
2123 SCMDCHECK(i!=-7,"error allocating memory");
2124 line2=strnword(line2,2);
2125 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
2126 SCMDCHECK(er!=-1,"file or data name not recognized");
2127 invk=cmd?cmd->invoke:0;
2128
2129 if(i!=-4) {
2130 inscan=1;
2131 molscancmd(sim,i,index,ms,cmd,cmdlistmols3);
2132 inscan=0; }
2133
2134 scmdflush(fptr);
2135 return CMDok;
2136
2137 scanportion:
2138 mptr=(moleculeptr) line2;
2139 scmdfprintf(cmd->cmds,fptr,"%i%,%i%,%i",invk,mptr->ident,mptr->mstate);
2140 scmdappenddata(cmd->cmds,dataid,1,3,(double)invk,(double)(mptr->ident),(double)(mptr->mstate));
2141 for(d=0;d<sim->dim;d++) {
2142 scmdfprintf(cmd->cmds,fptr,"%,%g",mptr->pos[d]);
2143 scmdappenddata(cmd->cmds,dataid,0,1,mptr->pos[d]); }
2144 scmdfprintf(cmd->cmds,fptr,"%,%s\n",molserno2string(mptr->serno,string));
2145 scmdappenddata(cmd->cmds,dataid,0,1,(double)(mptr->serno));
2146 return CMDok; }
2147
2148
2149 /* cmdlistmols4 */
cmdlistmols4(simptr sim,cmdptr cmd,char * line2)2150 enum CMDcode cmdlistmols4(simptr sim,cmdptr cmd,char *line2) {
2151 int i,d,*index,er;
2152 moleculeptr mptr;
2153 enum MolecState ms;
2154 static FILE *fptr;
2155 static int inscan=0,invk,dataid=-1;
2156 char string[STRCHAR];
2157
2158 if(inscan) goto scanportion;
2159 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
2160
2161 i=molstring2index1(sim,line2,&ms,&index);
2162 SCMDCHECK(i!=-1,"species is missing or cannot be read");
2163 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
2164 SCMDCHECK(i!=-3,"cannot read molecule state value");
2165 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
2166 SCMDCHECK(i!=-7,"error allocating memory");
2167 line2=strnword(line2,2);
2168 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
2169 SCMDCHECK(er!=-1,"file or data name not recognized");
2170 invk=cmd?cmd->invoke:0;
2171
2172 if(i!=-4) {
2173 inscan=1;
2174 molscancmd(sim,i,index,ms,cmd,cmdlistmols4);
2175 inscan=0; }
2176
2177 scmdflush(fptr);
2178 return CMDok;
2179
2180 scanportion:
2181 mptr=(moleculeptr) line2;
2182 scmdfprintf(cmd->cmds,fptr,"%i%,%i%,%i",invk,mptr->ident,mptr->mstate);
2183 scmdappenddata(cmd->cmds,dataid,1,3,(double)invk,(double)(mptr->ident),(double)(mptr->mstate));
2184 for(d=0;d<sim->dim;d++) {
2185 scmdfprintf(cmd->cmds,fptr,"%,%g",mptr->pos[d]+mptr->posoffset[d]);
2186 scmdappenddata(cmd->cmds,dataid,0,1,mptr->pos[d]+mptr->posoffset[d]); }
2187 scmdfprintf(cmd->cmds,fptr,"%,%s\n",molserno2string(mptr->serno,string));
2188 scmdappenddata(cmd->cmds,dataid,0,1,(double)(mptr->serno));
2189 return CMDok; }
2190
2191
2192 /* cmdlistmolscmpt */
cmdlistmolscmpt(simptr sim,cmdptr cmd,char * line2)2193 enum CMDcode cmdlistmolscmpt(simptr sim,cmdptr cmd,char *line2) {
2194 int i,*index,c,itct,d,er;
2195 moleculeptr mptr;
2196 enum MolecState ms;
2197 char cname[STRCHAR],string[STRCHAR];
2198 compartssptr cmptss;
2199 static FILE *fptr;
2200 static compartptr cmpt;
2201 static int inscan=0,invk,dataid=-1;
2202
2203 if(inscan) goto scanportion;
2204 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
2205
2206 i=molstring2index1(sim,line2,&ms,&index);
2207 SCMDCHECK(i!=-1,"species is missing or cannot be read");
2208 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
2209 SCMDCHECK(i!=-3,"cannot read molecule state value");
2210 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
2211 SCMDCHECK(i!=-7,"error allocating memory");
2212 line2=strnword(line2,2);
2213 SCMDCHECK(line2,"missing compartment name");
2214 itct=sscanf(line2,"%s",cname);
2215 SCMDCHECK(itct==1,"cannot read compartment name");
2216 cmptss=sim->cmptss;
2217 SCMDCHECK(cmptss,"no compartments defined");
2218 c=stringfind(cmptss->cnames,cmptss->ncmpt,cname);
2219 SCMDCHECK(c>=0,"compartment name not recognized");
2220 cmpt=cmptss->cmptlist[c];
2221 line2=strnword(line2,2);
2222 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
2223 SCMDCHECK(er!=-1,"file or data name not recognized");
2224 invk=cmd?cmd->invoke:0;
2225
2226 if(i!=-4) {
2227 inscan=1;
2228 molscancmd(sim,i,index,ms,cmd,cmdlistmolscmpt);
2229 inscan=0; }
2230
2231 scmdflush(fptr);
2232 return CMDok;
2233
2234 scanportion:
2235 mptr=(moleculeptr) line2;
2236 if(posincompart(sim,mptr->pos,cmpt,0)) {
2237 scmdfprintf(cmd->cmds,fptr,"%i%,%i%,%i",invk,mptr->ident,mptr->mstate);
2238 scmdappenddata(cmd->cmds,dataid,1,3,(double)invk,(double)(mptr->ident),(double)(mptr->mstate));
2239 for(d=0;d<sim->dim;d++) {
2240 scmdfprintf(cmd->cmds,fptr,"%,%g",mptr->pos[d]);
2241 scmdappenddata(cmd->cmds,dataid,0,1,mptr->pos[d]); }
2242 scmdfprintf(cmd->cmds,fptr,"%,%s\n",molserno2string(mptr->serno,string));
2243 scmdappenddata(cmd->cmds,dataid,0,1,(double)(mptr->serno)); }
2244 return CMDok; }
2245
2246
2247 /* cmdmolpos */
cmdmolpos(simptr sim,cmdptr cmd,char * line2)2248 enum CMDcode cmdmolpos(simptr sim,cmdptr cmd,char *line2) {
2249 int i,d,*index,er;
2250 moleculeptr mptr;
2251 enum MolecState ms;
2252 static FILE *fptr;
2253 static int inscan=0,dataid=-1;
2254
2255 if(inscan) goto scanportion;
2256 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
2257
2258 i=molstring2index1(sim,line2,&ms,&index);
2259 SCMDCHECK(i!=-1,"species is missing or cannot be read");
2260 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
2261 SCMDCHECK(i!=-3,"cannot read molecule state value");
2262 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
2263 SCMDCHECK(i!=-7,"error allocating memory");
2264 line2=strnword(line2,2);
2265 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
2266 SCMDCHECK(er!=-1,"file or data name not recognized");
2267
2268 scmdfprintf(cmd->cmds,fptr,"%g",sim->time);
2269 scmdappenddata(cmd->cmds,dataid,1,1,sim->time);
2270 if(i!=-4) {
2271 inscan=1;
2272 molscancmd(sim,i,index,ms,cmd,cmdmolpos);
2273 inscan=0; }
2274
2275 scmdfprintf(cmd->cmds,fptr,"\n");
2276 scmdflush(fptr);
2277 return CMDok;
2278
2279 scanportion:
2280 mptr=(moleculeptr) line2;
2281 for(d=0;d<sim->dim;d++) {
2282 scmdfprintf(cmd->cmds,fptr,"%,%g",mptr->pos[d]);
2283 scmdappenddata(cmd->cmds,dataid,0,1,mptr->pos[d]); }
2284 return CMDok; }
2285
2286
2287 /* cmdtrackmol */
cmdtrackmol(simptr sim,cmdptr cmd,char * line2)2288 enum CMDcode cmdtrackmol(simptr sim,cmdptr cmd,char *line2) {
2289 int itct,d,c,er;
2290 moleculeptr mptr;
2291 char string[STRCHAR];
2292 static FILE *fptr;
2293 static unsigned long long serno;
2294 static int inscan=0,dataid=-1;
2295
2296 if(inscan) goto scanportion;
2297 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
2298
2299 itct=sscanf(line2,"%s",string);
2300 SCMDCHECK(itct==1,"cannot read molecule serial number");
2301 serno=molstring2serno(string);
2302 SCMDCHECK(serno>0,"cannot read molecule serial number");
2303 line2=strnword(line2,2);
2304 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
2305 SCMDCHECK(er!=-1,"file or data name not recognized");
2306
2307 inscan=1;
2308 molscancmd(sim,-1,NULL,MSall,cmd,cmdtrackmol);
2309 inscan=0;
2310
2311 scmdflush(fptr);
2312 return CMDok;
2313
2314 scanportion:
2315 mptr=(moleculeptr) line2;
2316 if(!(mptr->serno==serno || (serno<0xFFFFFFFF && (mptr->serno&0xFFFFFFFF)==serno) || (serno<0xFFFFFFFF && mptr->serno>0xFFFFFFFF && (mptr->serno)>>32==serno)))
2317 return CMDok;
2318 scmdfprintf(cmd->cmds,fptr,"%g%,%s%,%s",sim->time,sim->mols->spname[mptr->ident],molms2string(mptr->mstate,string));
2319 scmdappenddata(cmd->cmds,dataid,1,3,sim->time,(double)(mptr->ident),(double)(mptr->mstate));
2320 scmdfprintf(cmd->cmds,fptr,"%,%s",molserno2string(mptr->serno,string));
2321 scmdappenddata(cmd->cmds,dataid,0,1,(double)(mptr->serno));
2322 for(d=0;d<sim->dim;d++) {
2323 scmdfprintf(cmd->cmds,fptr,"%,%g",mptr->pos[d]);
2324 scmdappenddata(cmd->cmds,dataid,0,1,mptr->pos[d]); }
2325 if(sim->cmptss)
2326 for(c=0;c<sim->cmptss->ncmpt;c++) {
2327 if(posincompart(sim,mptr->pos,sim->cmptss->cmptlist[c],0)) {
2328 scmdfprintf(cmd->cmds,fptr,"%,in");
2329 scmdappenddata(cmd->cmds,dataid,0,1,(double)1.0); }
2330 else {
2331 scmdfprintf(cmd->cmds,fptr,"%,out");
2332 scmdappenddata(cmd->cmds,dataid,0,1,(double)0.0); }}
2333 scmdfprintf(cmd->cmds,fptr,"\n");
2334 return CMDstop; }
2335
2336
2337 /* cmdmolmoments */
cmdmolmoments(simptr sim,cmdptr cmd,char * line2)2338 enum CMDcode cmdmolmoments(simptr sim,cmdptr cmd,char *line2) {
2339 int i,*index,d,d2,dim,er,dataid;
2340 FILE *fptr;
2341 moleculeptr mptr;
2342 enum MolecState ms;
2343 static double v1[DIMMAX],m1[DIMMAX*DIMMAX];
2344 static int inscan=0,ctr;
2345
2346 if(inscan) goto scanportion;
2347 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
2348
2349 i=molstring2index1(sim,line2,&ms,&index);
2350 SCMDCHECK(i!=-1,"species is missing or cannot be read");
2351 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
2352 SCMDCHECK(i!=-3,"cannot read molecule state value");
2353 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
2354 SCMDCHECK(i!=-7,"error allocating memory");
2355 line2=strnword(line2,2);
2356 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
2357 SCMDCHECK(er!=-1,"file or data name not recognized");
2358
2359 dim=sim->dim;
2360 ctr=0;
2361 for(d=0;d<dim;d++) v1[d]=0;
2362 for(d=0;d<dim*dim;d++) m1[d]=0;
2363
2364 if(i!=-4) {
2365 inscan=1;
2366 molscancmd(sim,i,index,ms,cmd,cmdmolmoments);
2367 for(d=0;d<dim;d++) v1[d]/=1.0*ctr;
2368 inscan=2;
2369 molscancmd(sim,i,index,ms,cmd,cmdmolmoments);
2370 inscan=0; }
2371
2372 scmdfprintf(cmd->cmds,fptr,"%g%,%i",sim->time,ctr);
2373 scmdappenddata(cmd->cmds,dataid,1,2,sim->time,(double)ctr);
2374 for(d=0;d<dim;d++) {
2375 scmdfprintf(cmd->cmds,fptr,"%,%g",v1[d]);
2376 scmdappenddata(cmd->cmds,dataid,0,1,v1[d]); }
2377 for(d=0;d<dim;d++)
2378 for(d2=0;d2<dim;d2++) {
2379 scmdfprintf(cmd->cmds,fptr,"%,%g",m1[d*dim+d2]/ctr);
2380 scmdappenddata(cmd->cmds,dataid,0,1,m1[d*dim+d2]/ctr); }
2381 scmdfprintf(cmd->cmds,fptr,"\n");
2382 scmdflush(fptr);
2383 return CMDok;
2384
2385 scanportion:
2386 mptr=(moleculeptr) line2;
2387 if(inscan==1) {
2388 ctr++;
2389 for(d=0;d<sim->dim;d++) v1[d]+=mptr->pos[d]; }
2390 else {
2391 for(d=0;d<sim->dim;d++)
2392 for(d2=0;d2<sim->dim;d2++)
2393 m1[d*sim->dim+d2]+=(mptr->pos[d]-v1[d])*(mptr->pos[d2]-v1[d2]); }
2394 return CMDok; }
2395
2396
2397 /* cmdsavesim */
cmdsavesim(simptr sim,cmdptr cmd,char * line2)2398 enum CMDcode cmdsavesim(simptr sim,cmdptr cmd,char *line2) {
2399 int er;
2400 FILE *fptr;
2401
2402 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
2403 er=scmdgetfptr(sim->cmds,line2,1,&fptr,NULL);
2404 SCMDCHECK(er!=-1,"file name not recognized");
2405 if(line2) {
2406 strcutwhite(line2,2); }
2407
2408 scmdfprintf(cmd->cmds,fptr,"# Configuration file automatically created by Smoldyn\n\n");
2409 writesim(sim,fptr);
2410 writegraphss(sim,fptr);
2411 writemols(sim,fptr);
2412 writewalls(sim,fptr);
2413 writesurfaces(sim,fptr);
2414 writecomparts(sim,fptr);
2415 writereactions(sim,fptr);
2416 writerules(sim,fptr);
2417 writelattices(sim,fptr);
2418 scmdwritecommands(sim->cmds,fptr,line2);
2419 writemolecules(sim,fptr);
2420 scmdfprintf(cmd->cmds,fptr,"\nend_file\n");
2421 scmdflush(fptr);
2422 return CMDok; }
2423
2424
2425 /* cmdmeansqrdispfree */
cmdmeansqrdispfree(cmdptr cmd)2426 void cmdmeansqrdispfree(cmdptr cmd) {
2427 int j;
2428
2429 if(cmd->v2 && cmd->i1)
2430 for(j=0;j<cmd->i1;j++)
2431 free(((double**)(cmd->v2))[j]);
2432 free(cmd->v2);
2433 free(cmd->v1);
2434 return; }
2435
2436
2437 /* cmdmeansqrdisp */
cmdmeansqrdisp(simptr sim,cmdptr cmd,char * line2)2438 enum CMDcode cmdmeansqrdisp(simptr sim,cmdptr cmd,char *line2) {
2439 char dimstr[STRCHAR];
2440 int i,*index,j,d,itct,dim,er,dataid;
2441 FILE *fptr;
2442 double r2,diff,**v2;
2443 long int *v1;
2444 enum MolecState ms;
2445 moleculeptr mptr;
2446 static double sum,sum4;
2447 static int inscan=0,ctr,msddim;
2448
2449 if(inscan) goto scanportion;
2450 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
2451
2452 i=molstring2index1(sim,line2,&ms,&index);
2453 SCMDCHECK(i!=-1,"species is missing or cannot be read");
2454 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
2455 SCMDCHECK(i!=-3,"cannot read molecule state value");
2456 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
2457 SCMDCHECK(i!=-7,"error allocating memory");
2458 line2=strnword(line2,2);
2459 SCMDCHECK(line2,"missing dimension information");
2460 itct=sscanf(line2,"%s",dimstr);
2461 SCMDCHECK(itct==1,"cannot read dimension information");
2462 msddim=1;
2463 if(!strcmp(dimstr,"all")) msddim=-1;
2464 else if(!strcmp(dimstr,"0") || !strcmp(dimstr,"x")) msddim=0;
2465 else if(!strcmp(dimstr,"1") || !strcmp(dimstr,"y")) msddim=1;
2466 else if(!strcmp(dimstr,"2") || !strcmp(dimstr,"z")) msddim=2;
2467 else msddim=3;
2468 SCMDCHECK(msddim<sim->dim,"invalid dimension value");
2469 line2=strnword(line2,2);
2470 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
2471 SCMDCHECK(er!=-1,"file or data name not recognized");
2472
2473 SCMDCHECK(cmd->i2!=2,"error on setup"); // failed before, don't try again
2474
2475 dim=sim->dim;
2476
2477 if(!cmd->i2) { // test for first run
2478 cmd->i2=1; // if first run, set up data structures
2479 ctr=(i==-4)?0:molcount(sim,i,index,ms,-1);
2480 cmd->i1=ctr; // size of arrays
2481 SCMDCHECK(ctr>0,"no molecules to track");
2482 cmd->freefn=&cmdmeansqrdispfree;
2483 cmd->v1=calloc(ctr,sizeof(long int)); // v1 is serial numbers
2484 if(!cmd->v1) {cmd->i2=2;return CMDwarn;}
2485 for(j=0;j<ctr;j++) ((long int*)cmd->v1)[j]=0;
2486 cmd->v2=calloc(ctr,sizeof(double**)); // v2 is positions
2487 if(!cmd->v2) {cmd->i2=2;return CMDwarn;}
2488 for(j=0;j<ctr;j++) ((double**)cmd->v2)[j]=NULL;
2489 for(j=0;j<ctr;j++) {
2490 ((double**)cmd->v2)[j]=(double*)calloc(dim,sizeof(double));
2491 if(!((double**)cmd->v2)[j]) {cmd->i2=2;return CMDwarn;}
2492 for(d=0;d<dim;d++) ((double**)cmd->v2)[j][d]=0; }
2493 ctr=0;
2494 if(i!=-4) {
2495 inscan=1;
2496 molscancmd(sim,i,index,ms,cmd,cmdmeansqrdisp);
2497 inscan=0; }
2498 sortVliv((long int*)cmd->v1,(void**)cmd->v2,cmd->i1); }
2499
2500 ctr=0; // start of code that is run every invocation
2501 sum=0;
2502 sum4=0;
2503
2504 if(i!=-4) {
2505 inscan=2;
2506 molscancmd(sim,i,index,ms,cmd,cmdmeansqrdisp);
2507 inscan=0; }
2508
2509 scmdfprintf(cmd->cmds,fptr,"%g%,%g%,%g\n",sim->time,sum/ctr,sum4/ctr);
2510 scmdappenddata(cmd->cmds,dataid,1,3,sim->time,sum/ctr,sum4/ctr);
2511
2512 scmdflush(fptr);
2513 return CMDok;
2514
2515 scanportion:
2516 mptr=(moleculeptr) line2;
2517 if(inscan==1) {
2518 ((long int*)cmd->v1)[ctr]=(long int)(mptr->serno&0xFFFFFFFF);
2519 for(d=0;d<sim->dim;d++)
2520 ((double**)cmd->v2)[ctr][d]=mptr->posoffset[d]+mptr->pos[d];
2521 ctr++; }
2522 else {
2523 v1=(long int*)cmd->v1;
2524 v2=(double**)cmd->v2;
2525 j=locateVli(v1,(long int)(mptr->serno&0xFFFFFFFF),cmd->i1);
2526 if(j>=0) {
2527 ctr++;
2528 if(msddim<0) {
2529 r2=0;
2530 for(d=0;d<sim->dim;d++) {
2531 diff=mptr->posoffset[d]+mptr->pos[d]-v2[j][d];
2532 r2+=diff*diff; }
2533 sum+=r2;
2534 sum4+=r2*r2; }
2535 else {
2536 diff=mptr->posoffset[msddim]+mptr->pos[msddim]-v2[j][msddim];
2537 sum+=diff*diff;
2538 sum4+=diff*diff*diff*diff; }}}
2539 return CMDok; }
2540
2541
2542 /* cmdmeansqrdisp2 */
cmdmeansqrdisp2(simptr sim,cmdptr cmd,char * line2)2543 enum CMDcode cmdmeansqrdisp2(simptr sim,cmdptr cmd,char *line2) {
2544 char dimstr[STRCHAR];
2545 int i,j,d,itct,dim,msddim,maxmoment,mom,*index,er,dataid;
2546 FILE *fptr;
2547 moleculeptr mptr;
2548 double sum[17];
2549 double r2,diff,**v2,*dblptr;
2550 long int *v1;
2551 enum MolecState ms;
2552 char reportchar;
2553 static char startchar;
2554 static int inscan=0,maxmol,ctr;
2555
2556 if(inscan) goto scanportion;
2557 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
2558
2559 i=molstring2index1(sim,line2,&ms,&index);
2560 SCMDCHECK(i!=-1,"species is missing or cannot be read");
2561 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
2562 SCMDCHECK(i!=-3,"cannot read molecule state value");
2563 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
2564 SCMDCHECK(i!=-7,"error allocating memory");
2565 line2=strnword(line2,2);
2566 SCMDCHECK(line2,"missing dimension information");
2567 itct=sscanf(line2,"%s",dimstr);
2568 SCMDCHECK(itct==1,"cannot read dimension information");
2569 if(!strcmp(dimstr,"all")) msddim=-1;
2570 else if(!strcmp(dimstr,"0") || !strcmp(dimstr,"x")) msddim=0;
2571 else if(!strcmp(dimstr,"1") || !strcmp(dimstr,"y")) msddim=1;
2572 else if(!strcmp(dimstr,"2") || !strcmp(dimstr,"z")) msddim=2;
2573 else msddim=3;
2574 SCMDCHECK(msddim<sim->dim,"invalid dimension value");
2575 line2=strnword(line2,2);
2576 SCMDCHECK(line2,"insufficient arguments");
2577 itct=strmathsscanf(line2,"%c %c %mi %mi",Varnames,Varvalues,Nvar,&startchar,&reportchar,&maxmol,&maxmoment);
2578 SCMDCHECK(itct==4,"cannot read start, report, max_mol, or max_moment information");
2579 SCMDCHECK(maxmol>0,"max_mol has to be at least 1");
2580 SCMDCHECK(maxmoment>0,"maxmoment has to be at least 1");
2581 SCMDCHECK(maxmoment<=16,"max_moment cannot exceed 16");
2582 line2=strnword(line2,5);
2583 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
2584 SCMDCHECK(er!=-1,"file or data name not recognized");
2585
2586 SCMDCHECK(cmd->i2!=2,"error on setup"); // failed before, don't try again
2587
2588 dim=sim->dim;
2589
2590 if(!cmd->i2) { // test for first run
2591 cmd->i2=1; // if first run, set up data structures
2592 cmd->i1=maxmol;
2593 cmd->i3=0;
2594 cmd->freefn=&cmdmeansqrdispfree;
2595 cmd->v1=calloc(maxmol,sizeof(long int)); // v1 is serial numbers
2596 if(!cmd->v1) {cmd->i2=2;return CMDwarn;}
2597 v1=(long int*)cmd->v1;
2598 for(j=0;j<maxmol;j++) v1[j]=0;
2599 cmd->v2=calloc(maxmol,sizeof(double**)); // v2 is positions
2600 if(!cmd->v2) {cmd->i2=2;return CMDwarn;}
2601 v2=(double**)cmd->v2;
2602 for(j=0;j<maxmol;j++) v2[j]=NULL;
2603 for(j=0;j<maxmol;j++) {
2604 v2[j]=(double*)calloc(2*dim+1,sizeof(double));
2605 if(!v2[j]) {cmd->i2=2;return CMDwarn;}
2606 for(d=0;d<2*dim+1;d++) v2[j][d]=0; }
2607 ctr=0;
2608 if(i!=-4) {
2609 inscan=1;
2610 molscancmd(sim,i,index,ms,cmd,cmdmeansqrdisp2);
2611 inscan=0; }
2612 SCMDCHECK(ctr<maxmol,"insufficient allocated space");
2613 cmd->i3=ctr;
2614 if(cmd->i3>0) sortVliv(v1,(void**)cmd->v2,cmd->i3); }
2615
2616 v1=(long int*)cmd->v1; // start of code that is run every invocation
2617 v2=(double**)cmd->v2;
2618
2619 if(i!=-4) {
2620 inscan=2; // update tracking information for all tracked molecules
2621 molscancmd(sim,i,index,ms,cmd,cmdmeansqrdisp2);
2622 inscan=0; }
2623 if(cmd->i3==cmd->i1) {SCMDCHECK(0,"not enough allocated space");}
2624 if(startchar!='i') // resort lists if appropriate
2625 if(cmd->i3>0) sortVliv(v1,(void**)cmd->v2,cmd->i3);
2626
2627 for(mom=0;mom<=maxmoment;mom++)
2628 sum[mom]=0;
2629 for(j=0;j<cmd->i3;j++) { // find moments of all reported results
2630 if((reportchar=='e' && v2[j][0]==3.0) || (reportchar=='r' && v2[j][0]==2.0)) { // molecule should be recorded
2631 if(msddim<0) {
2632 r2=0;
2633 for(d=0;d<dim;d++) {
2634 diff=v2[j][dim+1+d]-v2[j][1+d];
2635 r2+=diff*diff; }
2636 r2=sqrt(r2);
2637 for(mom=0;mom<=maxmoment;mom++)
2638 sum[mom]+=pow(r2,mom); }
2639 else {
2640 diff=v2[j][dim+1+msddim]-v2[j][1+msddim];
2641 for(mom=0;mom<=maxmoment;mom++)
2642 sum[mom]+=pow(diff,mom); }}}
2643
2644 if(sum[0]>0) {
2645 scmdfprintf(cmd->cmds,fptr,"%g%,%g",sim->time,sum[0]); // display results
2646 scmdappenddata(cmd->cmds,dataid,1,2,sim->time,sum[0]);
2647 for(mom=1;mom<=maxmoment;mom++) {
2648 scmdfprintf(cmd->cmds,fptr,"%,%g",sum[mom]/sum[0]);
2649 scmdappenddata(cmd->cmds,dataid,0,1,sum[mom]/sum[0]); }
2650 scmdfprintf(cmd->cmds,fptr,"\n"); }
2651
2652 for(j=0;j<cmd->i3;j++) { // stop tracking expired molecules
2653 if(v2[j][0]==0 || v2[j][0]==2.0) {
2654 v1[j]=v1[cmd->i3-1];
2655 v1[cmd->i3-1]=0;
2656 dblptr=v2[j];
2657 v2[j]=v2[cmd->i3-1];
2658 v2[cmd->i3-1]=dblptr;
2659 v2[cmd->i3-1][0]=0;
2660 j--;
2661 cmd->i3--; }
2662 else
2663 v2[j][0]-=1.0; }
2664 if(cmd->i3>0) sortVliv(v1,(void**)cmd->v2,cmd->i3);
2665
2666 scmdflush(fptr);
2667 return CMDok;
2668
2669 scanportion:
2670 mptr=(moleculeptr) line2;
2671 v1=(long int*)cmd->v1;
2672 v2=(double**)cmd->v2;
2673 if(inscan==1) {
2674 if(ctr==maxmol) return CMDstop;
2675 v1[ctr]=(long int)(mptr->serno&0xFFFFFFFF);
2676 if(startchar=='c') v2[ctr][0]=0;
2677 else v2[ctr][0]=2.0;
2678 for(d=0;d<sim->dim;d++)
2679 v2[ctr][1+d]=v2[ctr][sim->dim+1+d]=mptr->posoffset[d]+mptr->pos[d];
2680 ctr++; }
2681 else {
2682 j=locateVli(v1,(long int)(mptr->serno&0xFFFFFFFF),cmd->i3);
2683 if(j>=0) { // molecule was found
2684 v2[j][0]+=1.0;
2685 if(v2[j][0]==3.0) { // molecule is being tracked and exists, so record current positions
2686 for(d=0;d<sim->dim;d++)
2687 v2[j][sim->dim+1+d]=mptr->posoffset[d]+mptr->pos[d]; }}
2688 else if(startchar!='i') { // molecule was not found but should be tracked
2689 if(cmd->i3==cmd->i1) return CMDstop;
2690 j=cmd->i3++; // find empty spot
2691 v1[j]=(long int)(mptr->serno&0xFFFFFFFF);
2692 v2[j][0]=3.0;
2693 for(d=0;d<sim->dim;d++)
2694 v2[j][1+d]=v2[j][sim->dim+1+d]=mptr->posoffset[d]+mptr->pos[d]; }}
2695 return CMDok; }
2696
2697
2698 /* cmdmeansqrdisp3 */
cmdmeansqrdisp3(simptr sim,cmdptr cmd,char * line2)2699 enum CMDcode cmdmeansqrdisp3(simptr sim,cmdptr cmd,char *line2) {
2700 char dimstr[STRCHAR];
2701 int i,j,d,itct,dim,msddim,*index,er,dataid;
2702 FILE *fptr;
2703 moleculeptr mptr;
2704 double sum,wgt;
2705 double r2,diff,**v2,*dblptr,change;
2706 long int *v1;
2707 enum MolecState ms;
2708 char reportchar;
2709 static char startchar;
2710 static int inscan=0,ctr,maxmol;
2711
2712 if(inscan) goto scanportion;
2713 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
2714
2715 i=molstring2index1(sim,line2,&ms,&index);
2716 SCMDCHECK(i!=-1,"species is missing or cannot be read");
2717 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
2718 SCMDCHECK(i!=-3,"cannot read molecule state value");
2719 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
2720 SCMDCHECK(i!=-7,"error allocating memory");
2721 line2=strnword(line2,2);
2722 SCMDCHECK(line2,"missing dimension information");
2723 itct=sscanf(line2,"%s",dimstr);
2724 SCMDCHECK(itct==1,"cannot read dimension information");
2725 if(!strcmp(dimstr,"all")) msddim=-1;
2726 else if(!strcmp(dimstr,"0") || !strcmp(dimstr,"x")) msddim=0;
2727 else if(!strcmp(dimstr,"1") || !strcmp(dimstr,"y")) msddim=1;
2728 else if(!strcmp(dimstr,"2") || !strcmp(dimstr,"z")) msddim=2;
2729 else msddim=3;
2730 SCMDCHECK(msddim<sim->dim,"invalid dimension value");
2731 line2=strnword(line2,2);
2732 SCMDCHECK(line2,"insufficient arguments");
2733 itct=strmathsscanf(line2,"%c %c %i %mlg",Varnames,Varvalues,Nvar,&startchar,&reportchar,&maxmol,&change);
2734 SCMDCHECK(itct==4,"cannot read start, report, max_mol, or change information");
2735 SCMDCHECK(maxmol>0,"max_mol has to be at least 1");
2736 line2=strnword(line2,5);
2737 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
2738 SCMDCHECK(er!=-1,"file or data name not recognized");
2739 line2=strnword(line2,2);
2740 SCMDCHECK(change<=0 || line2,"missing task to be accomplished if change is small");
2741
2742 SCMDCHECK(cmd->i2!=2,"error on setup"); // failed before, don't try again
2743
2744 dim=sim->dim;
2745
2746 if(!cmd->i2) { // test for first run
2747 cmd->i2=1; // if first run, set up data structures
2748 cmd->i1=maxmol;
2749 cmd->i3=0;
2750 cmd->f1=-1;
2751 cmd->freefn=&cmdmeansqrdispfree;
2752 cmd->v1=calloc(maxmol,sizeof(long int)); // v1 is serial numbers
2753 if(!cmd->v1) {cmd->i2=2;return CMDwarn;}
2754 v1=(long int*)cmd->v1;
2755 for(j=0;j<maxmol;j++) v1[j]=0;
2756 cmd->v2=calloc(maxmol,sizeof(double**)); // v2 is positions
2757 if(!cmd->v2) {cmd->i2=2;return CMDwarn;}
2758 v2=(double**)cmd->v2;
2759 for(j=0;j<maxmol;j++) v2[j]=NULL;
2760 for(j=0;j<maxmol;j++) {
2761 v2[j]=(double*)calloc(2*dim+2,sizeof(double));
2762 if(!v2[j]) {cmd->i2=2;return CMDwarn;}
2763 for(d=0;d<2*dim+2;d++) v2[j][d]=0; }
2764 ctr=0;
2765 if(i!=-4) {
2766 inscan=1;
2767 molscancmd(sim,i,index,ms,cmd,cmdmeansqrdisp3);
2768 inscan=0; }
2769 SCMDCHECK(ctr<maxmol,"insufficient allocated space");
2770 cmd->i3=ctr;
2771 if(cmd->i3>0) sortVliv(v1,(void**)cmd->v2,cmd->i3); }
2772
2773 v1=(long int*)cmd->v1; // start of code that is run every invocation
2774 v2=(double**)cmd->v2;
2775
2776 if(i!=-4) {
2777 inscan=2;
2778 molscancmd(sim,i,index,ms,cmd,cmdmeansqrdisp3);
2779 inscan=0; }
2780 if(cmd->i3==cmd->i1) {SCMDCHECK(0,"not enough allocated space");}
2781 if(startchar!='i') // resort lists if appropriate
2782 if(cmd->i3>0) sortVliv(v1,(void**)cmd->v2,cmd->i3);
2783
2784 sum=0;
2785 ctr=0;
2786 wgt=0;
2787 for(j=0;j<cmd->i3;j++) { // find effective diffusion coefficients of all reported results
2788 if((reportchar=='e' && v2[j][0]==3.0) || (reportchar=='r' && v2[j][0]==2.0)) { // molecule should be recorded
2789 ctr++;
2790 if(msddim<0) {
2791 r2=0;
2792 for(d=0;d<dim;d++) {
2793 diff=v2[j][dim+1+d]-v2[j][1+d];
2794 r2+=diff*diff; }
2795 sum+=r2;
2796 wgt+=sim->time-v2[j][2*dim+1]; }
2797 else {
2798 diff=v2[j][dim+1+msddim]-v2[j][1+msddim];
2799 r2=diff*diff;
2800 sum+=r2;
2801 wgt+=sim->time-v2[j][2*dim+1]; }}}
2802 if(msddim<0) sum/=(2.0*dim);
2803 else sum/=2.0;
2804
2805 scmdfprintf(cmd->cmds,fptr,"%g%,%i%,%g\n",sim->time,ctr,sum/wgt); // display results
2806 scmdappenddata(cmd->cmds,dataid,1,3,sim->time,(double)ctr,sum/wgt);
2807
2808 for(j=0;j<cmd->i3;j++) { // stop tracking expired molecules
2809 if(v2[j][0]==0 || v2[j][0]==2.0) {
2810 v1[j]=v1[cmd->i3-1];
2811 v1[cmd->i3-1]=0;
2812 dblptr=v2[j];
2813 v2[j]=v2[cmd->i3-1];
2814 v2[cmd->i3-1]=dblptr;
2815 v2[cmd->i3-1][0]=0;
2816 j--;
2817 cmd->i3--; }
2818 else
2819 v2[j][0]-=1.0; }
2820 if(cmd->i3>0) sortVliv(v1,(void**)cmd->v2,cmd->i3);
2821
2822 if(change>0 && ctr>0 && cmd->f1>0 && fabs((sum/ctr-cmd->f1)/cmd->f1)<change)
2823 return docommand(sim,cmd,line2);
2824 cmd->f1=sum/ctr;
2825 scmdflush(fptr);
2826 return CMDok;
2827
2828 scanportion:
2829 mptr=(moleculeptr) line2;
2830 v1=(long int*)cmd->v1;
2831 v2=(double**)cmd->v2;
2832 if(inscan==1) {
2833 if(ctr==maxmol) return CMDstop;
2834 v1[ctr]=(long int)(mptr->serno&0xFFFFFFFF);
2835 if(startchar=='c') v2[ctr][0]=0;
2836 else v2[ctr][0]=2.0;
2837 for(d=0;d<sim->dim;d++)
2838 v2[ctr][1+d]=v2[ctr][sim->dim+1+d]=mptr->posoffset[d]+mptr->pos[d];
2839 v2[ctr][2*sim->dim+1]=sim->time;
2840 ctr++; }
2841 else {
2842 j=locateVli(v1,(long int)(mptr->serno&0xFFFFFFFF),cmd->i3);
2843 if(j>=0) { // molecule was found
2844 v2[j][0]+=1.0;
2845 if(v2[j][0]==3.0) { // molecule is being tracked and exists, so record current positions
2846 for(d=0;d<sim->dim;d++)
2847 v2[j][sim->dim+1+d]=mptr->posoffset[d]+mptr->pos[d]; }}
2848 else if(startchar!='i') { // molecule was not found but should be tracked
2849 if(cmd->i3==cmd->i1) return CMDstop;
2850 j=cmd->i3++; // find empty spot
2851 v1[j]=(long int)(mptr->serno&0xFFFFFFFF);
2852 v2[j][0]=3.0;
2853 for(d=0;d<sim->dim;d++)
2854 v2[j][1+d]=v2[j][sim->dim+1+d]=mptr->posoffset[d]+mptr->pos[d];
2855 v2[j][2*sim->dim+1]=sim->time; }}
2856 return CMDok; }
2857
2858
2859 /* cmdresidencetime */
cmdresidencetime(simptr sim,cmdptr cmd,char * line2)2860 enum CMDcode cmdresidencetime(simptr sim,cmdptr cmd,char *line2) {
2861 int i,j,d,itct,summaryout,listout,*index,er,dataid;
2862 FILE *fptr;
2863 moleculeptr mptr;
2864 double sum;
2865 double **v2,*dblptr;
2866 long int *v1;
2867 enum MolecState ms;
2868 char reportchar;
2869 static char startchar;
2870 static int inscan=0,ctr,maxmol;
2871
2872 if(inscan) goto scanportion;
2873 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
2874
2875 i=molstring2index1(sim,line2,&ms,&index);
2876 SCMDCHECK(i!=-1,"species is missing or cannot be read");
2877 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
2878 SCMDCHECK(i!=-3,"cannot read molecule state value");
2879 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
2880 SCMDCHECK(i!=-7,"error allocating memory");
2881 line2=strnword(line2,2);
2882 SCMDCHECK(line2,"insufficient arguments");
2883 itct=strmathsscanf(line2,"%c %c %mi %mi %mi",Varnames,Varvalues,Nvar,&startchar,&reportchar,&summaryout,&listout,&maxmol);
2884 SCMDCHECK(itct==5,"cannot read start, report, summary_out, list_out, or max_mol information");
2885 SCMDCHECK(maxmol>0,"max_mol has to be at least 1");
2886 line2=strnword(line2,6);
2887 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
2888 SCMDCHECK(er!=-1,"file or data name not recognized");
2889
2890 SCMDCHECK(cmd->i2!=2,"error on setup"); // failed before, don't try again
2891
2892 if(!cmd->i2) { // test for first run
2893 cmd->i2=1; // if first run, set up data structures
2894 cmd->i1=maxmol;
2895 cmd->i3=0;
2896 cmd->freefn=&cmdmeansqrdispfree;
2897 cmd->v1=calloc(maxmol,sizeof(long int)); // v1 is serial numbers
2898 if(!cmd->v1) {cmd->i2=2;return CMDwarn;}
2899 v1=(long int*)cmd->v1;
2900 for(j=0;j<maxmol;j++) v1[j]=0;
2901 cmd->v2=calloc(maxmol,sizeof(double**)); // v2 is creation times
2902 if(!cmd->v2) {cmd->i2=2;return CMDwarn;}
2903 v2=(double**)cmd->v2;
2904 for(j=0;j<maxmol;j++) v2[j]=NULL;
2905 for(j=0;j<maxmol;j++) {
2906 v2[j]=(double*)calloc(2,sizeof(double));
2907 if(!v2[j]) {cmd->i2=2;return CMDwarn;}
2908 for(d=0;d<2;d++) v2[j][d]=0; }
2909 ctr=0;
2910 if(i!=-4) {
2911 inscan=1;
2912 molscancmd(sim,i,index,ms,cmd,cmdresidencetime);
2913 inscan=0; }
2914 SCMDCHECK(ctr<maxmol,"insufficient allocated space");
2915 cmd->i3=ctr;
2916 if(cmd->i3>0) sortVliv(v1,(void**)cmd->v2,cmd->i3); }
2917
2918 v1=(long int*)cmd->v1; // start of code that is run every invocation
2919 v2=(double**)cmd->v2;
2920
2921 if(i!=-4) {
2922 inscan=2;
2923 molscancmd(sim,i,index,ms,cmd,cmdresidencetime);
2924 inscan=0; }
2925
2926 if(cmd->i3==cmd->i1) {SCMDCHECK(0,"not enough allocated space");}
2927 if(startchar!='i') // resort lists if appropriate
2928 if(cmd->i3>0) sortVliv(v1,(void**)cmd->v2,cmd->i3);
2929
2930 sum=0;
2931 ctr=0;
2932 for(j=0;j<cmd->i3;j++) { // find effective diffusion coefficients of all reported results
2933 if((reportchar=='e' && v2[j][0]==3.0) || (reportchar=='r' && v2[j][0]==2.0)) { // molecule should be recorded
2934 ctr++;
2935 sum+=sim->time-v2[j][1];
2936 if(listout>0 && cmd->invoke>0 && cmd->invoke%listout==0) {
2937 scmdfprintf(cmd->cmds,fptr,"%li%,%g\n",v1[j],sim->time-v2[j][1]);
2938 scmdappenddata(cmd->cmds,dataid,1,2,(double)v1[j],sim->time-v2[j][1]); }}}
2939
2940 if(summaryout>0 && cmd->invoke>0 && cmd->invoke%summaryout==0) {
2941 scmdfprintf(cmd->cmds,fptr,"%g%,%i%,%g\n",sim->time,ctr,sum/(double)ctr); // display results
2942 scmdappenddata(cmd->cmds,dataid,1,3,sim->time,(double)ctr,sum/(double)ctr); }
2943
2944 for(j=0;j<cmd->i3;j++) { // stop tracking expired molecules
2945 if(v2[j][0]==0 || v2[j][0]==2.0) {
2946 v1[j]=v1[cmd->i3-1];
2947 v1[cmd->i3-1]=0;
2948 dblptr=v2[j];
2949 v2[j]=v2[cmd->i3-1];
2950 v2[cmd->i3-1]=dblptr;
2951 v2[cmd->i3-1][0]=0;
2952 j--;
2953 cmd->i3--; }
2954 else
2955 v2[j][0]-=1.0; }
2956 if(cmd->i3>0) sortVliv(v1,(void**)cmd->v2,cmd->i3);
2957
2958 scmdflush(fptr);
2959 return CMDok;
2960
2961 scanportion:
2962 mptr=(moleculeptr) line2;
2963 v1=(long int*)cmd->v1;
2964 v2=(double**)cmd->v2;
2965 if(inscan==1) {
2966 if(ctr==maxmol) return CMDstop;
2967 v1[ctr]=(long int)(mptr->serno&0xFFFFFFFF);
2968 if(startchar=='c') v2[ctr][0]=0;
2969 else v2[ctr][0]=2.0;
2970 v2[ctr][1]=sim->time;
2971 ctr++; }
2972 else {
2973 j=locateVli(v1,(long int)(mptr->serno&0xFFFFFFFF),cmd->i3);
2974 if(j>=0) { // molecule was found
2975 v2[j][0]+=1.0; }
2976 else if(startchar!='i') { // molecule was not found but should be tracked
2977 if(cmd->i3==cmd->i1) return CMDstop;
2978 j=cmd->i3++; // find empty spot
2979 v1[j]=(long int)(mptr->serno&0xFFFFFFFF);
2980 v2[j][0]=3.0;
2981 v2[j][1]=sim->time; }}
2982 return CMDok; }
2983
2984
2985 /* cmddiagnostics */
cmddiagnostics(simptr sim,cmdptr cmd,char * line2)2986 enum CMDcode cmddiagnostics(simptr sim,cmdptr cmd,char *line2) {
2987 int itct,order;
2988 static char nm[STRCHAR];
2989 enum SmolStruct ss;
2990
2991 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
2992 SCMDCHECK(line2,"missing argument");
2993 itct=sscanf(line2,"%s",nm);
2994 SCMDCHECK(itct==1,"read failure");
2995 ss=simstring2ss(nm);
2996 SCMDCHECK(ss!=SSnone,"diagnostic type not recognized");
2997
2998 if(ss==SSsim || ss==SSall) simoutput(sim);
2999 if(ss==SSwall || ss==SSall) walloutput(sim);
3000 if(ss==SSmolec || ss==SSall) molssoutput(sim);
3001 if(ss==SSsurf || ss==SSall) surfaceoutput(sim);
3002 if(ss==SScmd || ss==SSall) scmdoutput(sim->cmds);
3003 if(ss==SSbox || ss==SSall) boxssoutput(sim);
3004 if(ss==SSrxn || ss==SSall)
3005 for(order=0;order<MAXORDER;order++) rxnoutput(sim,order);
3006 if(ss==SSrule || ss==SSall) ruleoutput(sim);
3007 if(ss==SScmpt || ss==SSall) compartoutput(sim);
3008 if(ss==SSport || ss==SSall) portoutput(sim);
3009 if(ss==SScheck || ss==SSall) checksimparams(sim);
3010 return CMDok; }
3011
3012 /* cmdexecutiontime */
cmdexecutiontime(simptr sim,cmdptr cmd,char * line2)3013 enum CMDcode cmdexecutiontime(simptr sim,cmdptr cmd,char *line2) {
3014 int er,dataid;
3015 FILE *fptr;
3016
3017 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
3018 er=scmdgetfptr(sim->cmds,line2,3,&fptr,&dataid);
3019 SCMDCHECK(er!=-1,"file or data name not recognized");
3020
3021 scmdfprintf(cmd->cmds,fptr,"%g%,%g\n",sim->time,sim->elapsedtime+difftime(time(NULL),sim->clockstt));
3022 scmdappenddata(cmd->cmds,dataid,1,2,sim->time,sim->elapsedtime+difftime(time(NULL),sim->clockstt));
3023
3024 scmdflush(fptr);
3025 return CMDok; }
3026
3027
3028 /* cmdprintLattice */
cmdprintLattice(simptr sim,cmdptr cmd,char * line2)3029 enum CMDcode cmdprintLattice(simptr sim,cmdptr cmd,char *line2) {
3030 FILE *fptr;
3031 latticeptr lattice;
3032 char *buffer;
3033 int n,i,er;
3034
3035 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
3036 er=scmdgetfptr(sim->cmds,line2,1,&fptr,NULL);
3037 SCMDCHECK(er!=-1,"file name not recognized");
3038
3039 n=sim->latticess->nlattice;
3040 buffer=NULL;
3041 for(i=0;i<n;++i) {
3042 lattice=sim->latticess->latticelist[i];
3043 scmdfprintf(cmd->cmds,fptr,"Lattice %d: %s:\n",i,lattice->latticename);
3044 NSV_CALL(nsv_print(lattice->nsv,&buffer));
3045 scmdfprintf(cmd->cmds,fptr,"%s",buffer?buffer:"Error");
3046 buffer=NULL; }
3047 scmdflush(fptr);
3048 return CMDok; }
3049
3050
3051 /* cmdwriteVTK */
cmdwriteVTK(simptr sim,cmdptr cmd,char * line2)3052 enum CMDcode cmdwriteVTK(simptr sim,cmdptr cmd,char *line2) {
3053
3054 #ifndef OPTION_VTK
3055 (void)sim;
3056 (void)cmd;
3057 (void)line2;
3058 simLog(NULL,11,"ERROR: VTK option not set. Recompile with OPTION_VTK = ON\n");
3059 #else
3060
3061 char nm[STRCHAR];
3062 moleculeptr mptr;
3063 int itct;
3064 static vtkUnstructuredGrid* grid;
3065 static int inscan=0;
3066
3067 if(inscan) goto scanportion;
3068 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
3069
3070 SCMDCHECK(line2,"file name not given");
3071 itct=sscanf(line2,"%s",nm);
3072 SCMDCHECK(itct==1,"read failure");
3073
3074 // printf("Writing out VTK files....\n");
3075 grid = vtkCreateMolGrid();
3076
3077 inscan=1;
3078 molscancmd(sim,-1,NULL,MSall,cmd,cmdwriteVTK);
3079 inscan=0;
3080 vtkWriteGrid(nm,"Molecules",cmd->invoke,grid);
3081 vtkDeleteGrid(grid);
3082
3083 #ifdef OPTION_NSV
3084 static char nm2[STRCHAR];
3085 latticeptr lattice;
3086 int ll;
3087
3088 if (sim->latticess) {
3089 for (ll = 0; ll < sim->latticess->nlattice; ++ll) {
3090 lattice = sim->latticess->latticelist[ll];
3091 snprintf(nm2,STRCHAR,"Lattice%02d_",ll);
3092 switch(lattice->type) {
3093 case LATTICEnsv:
3094 if (lattice->nsv) vtkWriteGrid(nm,nm2,cmd->invoke,nsv_get_grid(lattice->nsv));
3095 break;
3096 case LATTICEpde:
3097 //not implemented
3098 break;
3099 case LATTICEnone:
3100 break; }}}
3101 #endif
3102
3103 return CMDok;
3104
3105 scanportion:
3106 mptr=(moleculeptr) line2;
3107 vtkAddPoint(grid,mptr->pos[0],mptr->pos[1],mptr->pos[2],(long int)(mptr->serno&0xFFFFFFFF),mptr->ident);
3108 #endif
3109 return CMDok; }
3110
3111
3112 /* cmdprintdata */
cmdprintdata(simptr sim,cmdptr cmd,char * line2)3113 enum CMDcode cmdprintdata(simptr sim,cmdptr cmd,char *line2) {
3114 FILE *fptr;
3115 int itct,er,did,dataid,i,j,erase;
3116 listptrdd list;
3117 cmdssptr cmds;
3118 char dname[STRCHAR];
3119
3120 if(line2 && !strcmp(line2,"cmdtype")) return CMDobserve;
3121
3122 cmds=sim->cmds;
3123 SCMDCHECK(line2,"missing data name");
3124 itct=sscanf(line2,"%s",dname);
3125 SCMDCHECK(itct==1,"cannot read data name");
3126 SCMDCHECK(cmds->ndata,"no data files have been declared");
3127 did=stringfind(cmds->dname,cmds->ndata,dname);
3128 SCMDCHECK(did>=0,"data name not recognized");
3129 list=cmds->data[did];
3130
3131 line2=strnword(line2,2);
3132 er=scmdgetfptr(cmds,line2,3,&fptr,&dataid);
3133 SCMDCHECK(er!=-1,"output file or data name not recognized");
3134
3135 erase=0;
3136 if(line2 && (line2=strnword(line2,2))) {
3137 itct=sscanf(line2,"%i",&erase);
3138 SCMDCHECK(itct==1,"erase value needs to be 0 or 1"); }
3139
3140 for(i=0;i<list->nrow;i++) {
3141 for(j=0;j<list->ncol;j++) {
3142 scmdfprintf(cmds,fptr,"%g",list->data[i*list->maxcol+j]);
3143 if(j<list->ncol-1) scmdfprintf(cmds,fptr,"%,");
3144 scmdappenddata(cmds,dataid,j==0?1:0,1,list->data[i*list->maxcol+j]); }
3145 scmdfprintf(cmds,fptr,"\n"); }
3146 scmdflush(fptr);
3147
3148 if(erase) ListClearDD(list);
3149
3150 return CMDok; }
3151
3152
3153 /**********************************************************/
3154 /****************** system manipulation ********************/
3155 /**********************************************************/
3156
3157
3158 /* cmdset */
cmdset(simptr sim,cmdptr cmd,char * line2)3159 enum CMDcode cmdset(simptr sim,cmdptr cmd,char *line2) {
3160 int er,itct;
3161 char word[STRCHAR];
3162
3163 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3164 SCMDCHECK(line2,"missing argument");
3165 itct=sscanf(line2,"%s",word);
3166 SCMDCHECK(itct==1,"missing statement");
3167 line2=strnword(line2,2);
3168 SCMDCHECK(line2,"missing statement text");
3169 er=simreadstring(sim,NULL,word,line2);
3170 SCMDCHECK(!er,"%s",ErrorString);
3171 return CMDok; }
3172
3173
3174 /* cmdpointsource */
cmdpointsource(simptr sim,cmdptr cmd,char * line2)3175 enum CMDcode cmdpointsource(simptr sim,cmdptr cmd,char *line2) {
3176 int itct,num,i;
3177 char nm[STRCHAR];
3178 double pos[DIMMAX];
3179
3180 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3181 SCMDCHECK(line2,"missing argument");
3182 SCMDCHECK(sim->mols,"molecules are undefined");
3183 itct=strmathsscanf(line2,"%s %mi",Varnames,Varvalues,Nvar,nm,&num);
3184 SCMDCHECK(itct==2,"read failure");
3185 SCMDCHECK(num>=0,"number cannot be negative");
3186 i=stringfind(sim->mols->spname,sim->mols->nspecies,nm);
3187 SCMDCHECK(i>=1,"name not recognized");
3188 line2=strnword(line2,3);
3189 SCMDCHECK(line2,"missing location");
3190 if(sim->dim==1) itct=strmathsscanf(line2,"%mlg",Varnames,Varvalues,Nvar,&pos[0]);
3191 else if(sim->dim==2) itct=strmathsscanf(line2,"%mlg %mlg",Varnames,Varvalues,Nvar,&pos[0],&pos[1]);
3192 else itct=strmathsscanf(line2,"%mlg %mlg %mlg",Varnames,Varvalues,Nvar,&pos[0],&pos[1],&pos[2]);
3193 SCMDCHECK(itct==sim->dim,"insufficient location dimensions");
3194 SCMDCHECK(addmol(sim,num,i,pos,pos,1)==0,"not enough available molecules");
3195 return CMDok; }
3196
3197
3198 /* cmdvolumesource */
cmdvolumesource(simptr sim,cmdptr cmd,char * line2)3199 enum CMDcode cmdvolumesource(simptr sim,cmdptr cmd,char *line2) {
3200 int itct,num,i,d;
3201 char nm[STRCHAR];
3202 double poslo[DIMMAX],poshi[DIMMAX],numdbl;
3203
3204 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3205 SCMDCHECK(line2,"missing argument");
3206 SCMDCHECK(sim->mols,"molecules are undefined");
3207 itct=strmathsscanf(line2,"%s %mlg",Varnames,Varvalues,Nvar,nm,&numdbl);
3208 SCMDCHECK(itct==2,"read failure");
3209 SCMDCHECK(numdbl>=0,"number cannot be negative");
3210 num=(int) numdbl;
3211 if(num!=numdbl) num=poisrandD(numdbl);
3212 i=stringfind(sim->mols->spname,sim->mols->nspecies,nm);
3213 SCMDCHECK(i>=1,"name not recognized");
3214 line2=strnword(line2,3);
3215 SCMDCHECK(line2,"missing location");
3216 for(d=0;d<sim->dim;d++) {
3217 SCMDCHECK(line2,"missing argument");
3218 itct=strmathsscanf(line2,"%mlg %mlg",Varnames,Varvalues,Nvar,&poslo[d],&poshi[d]);
3219 SCMDCHECK(itct==2,"read failure");
3220 line2=strnword(line2,3); }
3221 SCMDCHECK(addmol(sim,num,i,poslo,poshi,1)==0,"not enough available molecules");
3222 return CMDok; }
3223
3224
3225 /* guassiansource */
cmdgaussiansource(simptr sim,cmdptr cmd,char * line2)3226 enum CMDcode cmdgaussiansource(simptr sim,cmdptr cmd,char *line2) {
3227 int itct,num,i,d,imol,dim;
3228 char nm[STRCHAR];
3229 double mean[DIMMAX],sigma[DIMMAX],pos[DIMMAX],lowcorner[DIMMAX],highcorner[DIMMAX],numdbl;
3230
3231 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3232
3233 dim=sim->dim;
3234 SCMDCHECK(line2,"missing argument");
3235 SCMDCHECK(sim->mols,"molecules are undefined");
3236 itct=strmathsscanf(line2,"%s %mlg",Varnames,Varvalues,Nvar,nm,&numdbl);
3237 SCMDCHECK(itct==2,"read failure");
3238 SCMDCHECK(numdbl>=0,"number cannot be negative");
3239 num=(int) numdbl;
3240 if(num!=numdbl) num=poisrandD(numdbl);
3241 i=stringfind(sim->mols->spname,sim->mols->nspecies,nm);
3242 SCMDCHECK(i>=1,"name not recognized");
3243 line2=strnword(line2,3);
3244 SCMDCHECK(line2,"missing location");
3245 for(d=0;d<dim;d++) {
3246 SCMDCHECK(line2,"missing argument");
3247 itct=strmathsscanf(line2,"%mlg %mlg",Varnames,Varvalues,Nvar,&mean[d],&sigma[d]);
3248 SCMDCHECK(itct==2,"read failure");
3249 line2=strnword(line2,3); }
3250
3251 systemcorners(sim,lowcorner,highcorner);
3252 for(imol=0;imol<num;imol++) {
3253 for(d=0;d<dim;d++) {
3254 do {
3255 pos[d]=mean[d]+sigma[d]*gaussrandD(); }
3256 while(pos[d]<lowcorner[d] || pos[d]>highcorner[d]); }
3257 SCMDCHECK(addmol(sim,1,i,pos,pos,0)==0,"not enough available molecules"); }
3258 return CMDok; }
3259
3260
3261 /* cmdmovesurfacemol */
cmdmovesurfacemol(simptr sim,cmdptr cmd,char * line2)3262 enum CMDcode cmdmovesurfacemol(simptr sim,cmdptr cmd,char *line2) {
3263 int itct,i,s1,s2,d,*index;
3264 char nm[STRCHAR],nm2[STRCHAR];
3265 enum MolecState ms1;
3266 moleculeptr mptr;
3267 double pos[DIMMAX];
3268 static panelptr pnl2;
3269 static enum MolecState ms2;
3270 static surfaceptr srf,srf2;
3271 static double prob;
3272 static enum PanelShape ps1,ps2;
3273 static int inscan=0,p1,p2;
3274
3275 if(inscan) goto scanportion;
3276 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3277
3278 SCMDCHECK(line2,"missing arguments");
3279 SCMDCHECK(sim->mols,"molecules are undefined");
3280 SCMDCHECK(sim->srfss,"surfaces are undefined");
3281 itct=strmathsscanf(line2,"%s %mlg",Varnames,Varvalues,Nvar,nm,&prob);
3282 SCMDCHECK(itct==2,"failed to read molecule name or probability");
3283
3284 i=molstring2index1(sim,line2,&ms1,&index);
3285 SCMDCHECK(i!=-1,"species is missing or cannot be read");
3286 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
3287 SCMDCHECK(i!=-3,"cannot read molecule state value");
3288 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
3289 SCMDCHECK(i!=-7,"error allocating memory");
3290 SCMDCHECK(ms1==MSfront || ms1==MSback || ms1==MSup || ms1==MSdown || ms1==MSall,"illegal molecule state");
3291 SCMDCHECK(prob>=0 && prob<=1,"probability out of bounds");
3292 line2=strnword(line2,3);
3293 SCMDCHECK(line2,"missing originating surface:panel");
3294 itct=sscanf(line2,"%s %s",nm,nm2);
3295 SCMDCHECK(itct==2,"failed to read surfaces and panels");
3296 s1=readsurfacename(sim,nm,&ps1,&p1);
3297 SCMDCHECK(s1>=0,"failed to read surface1");
3298 SCMDCHECK(p1>=0 || p1==-5,"failed to read panel1");
3299 s2=readsurfacename(sim,nm2,&ps2,&p2);
3300 SCMDCHECK(s2>=0,"failed to read surface2");
3301 SCMDCHECK(p2>=0 || p2==-5,"failed to read panel2");
3302 line2=strnword(line2,3);
3303 if(line2) {
3304 itct=sscanf(line2,"%s",nm);
3305 SCMDCHECK(itct==1,"failed to read final state");
3306 ms2=molstring2ms(nm);
3307 SCMDCHECK(ms2!=MSnone,"failed to read final state");
3308 line2=strnword(line2,2); }
3309 else ms2=MSnone;
3310
3311 srf=sim->srfss->srflist[s1];
3312 srf2=sim->srfss->srflist[s2];
3313 if(p2==-5) pnl2=NULL;
3314 else pnl2=srf2->panels[ps2][p2];
3315
3316 if(i!=-4) {
3317 inscan=1;
3318 molscancmd(sim,i,index,ms1,cmd,cmdmovesurfacemol);
3319 inscan=0; }
3320
3321 sim->mols->touch++;
3322 return CMDok;
3323
3324 scanportion:
3325 mptr=(moleculeptr) line2;
3326 if(mptr->pnl && mptr->pnl->srf==srf && (p1==-5 || mptr->pnl==srf->panels[ps1][p1]))
3327 if(randCOD()<prob) {
3328 if(p2==-5) pnl2=surfrandpos(srf2,pos,sim->dim);
3329 else panelrandpos(pnl2,pos,sim->dim);
3330 for(d=0;d<sim->dim;d++) {
3331 mptr->posoffset[d]=mptr->pos[d]-pos[d];
3332 mptr->posx[d]=mptr->pos[d]=pos[d]; }
3333 molchangeident(sim,mptr,-1,-1,mptr->ident,ms2==MSnone?mptr->mstate:ms2,pnl2); }
3334 return CMDok; }
3335
3336
3337 /* cmdkillmol */
cmdkillmol(simptr sim,cmdptr cmd,char * line2)3338 enum CMDcode cmdkillmol(simptr sim,cmdptr cmd,char *line2) {
3339 int i,*index;
3340 moleculeptr mptr;
3341 enum MolecState ms;
3342 static int inscan=0;
3343
3344 if(inscan) goto scanportion;
3345 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3346
3347 i=molstring2index1(sim,line2,&ms,&index);
3348 SCMDCHECK(i!=-1,"species is missing or cannot be read");
3349 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
3350 SCMDCHECK(i!=-3,"cannot read molecule state value");
3351 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
3352 SCMDCHECK(i!=-7,"error allocating memory");
3353
3354 if(i!=-4) {
3355 inscan=1;
3356 molscancmd(sim,i,index,ms,cmd,cmdkillmol);
3357 inscan=0; }
3358 return CMDok;
3359
3360 scanportion:
3361 mptr=(moleculeptr) line2;
3362 molkill(sim,mptr,mptr->list,-1);
3363 return CMDok; }
3364
3365
3366 /* cmdkillmolprob */
cmdkillmolprob(simptr sim,cmdptr cmd,char * line2)3367 enum CMDcode cmdkillmolprob(simptr sim,cmdptr cmd,char *line2) {
3368 int itct,i,*index;
3369 moleculeptr mptr;
3370 enum MolecState ms;
3371 static double prob;
3372 static char probstr[STRCHAR];
3373 static int xyzvar;
3374 static int inscan=0;
3375
3376 if(inscan) goto scanportion;
3377 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3378
3379 i=molstring2index1(sim,line2,&ms,&index);
3380 SCMDCHECK(i!=-1,"species is missing or cannot be read");
3381 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
3382 SCMDCHECK(i!=-3,"cannot read molecule state value");
3383 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
3384 SCMDCHECK(i!=-7,"error allocating memory");
3385 line2=strnword(line2,2);
3386 SCMDCHECK(line2,"missing probability value");
3387 if(strhasname(line2,"x") || strhasname(line2,"y") || strhasname(line2,"z")) {
3388 xyzvar=1;
3389 strcpy(probstr,line2); }
3390 else {
3391 xyzvar=0;
3392 itct=strmathsscanf(line2,"%mlg",Varnames,Varvalues,Nvar,&prob);
3393 SCMDCHECK(itct==1,"killmolprob format: name[(state)] probability");
3394 SCMDCHECK(prob>=0 && prob<=1,"probability needs to be between 0 and 1"); }
3395
3396 if(i!=-4) {
3397 inscan=1;
3398 molscancmd(sim,i,index,ms,cmd,cmdkillmolprob);
3399 inscan=0; }
3400
3401 return CMDok;
3402
3403 scanportion:
3404 mptr=(moleculeptr) line2;
3405 if(xyzvar) {
3406 simsetvariable(sim,"x",mptr->pos[0]);
3407 if(sim->dim>1) simsetvariable(sim,"y",mptr->pos[1]);
3408 if(sim->dim>2) simsetvariable(sim,"z",mptr->pos[2]);
3409 strmathsscanf(probstr,"%mlg",Varnames,Varvalues,Nvar,&prob); }
3410 if(coinrandD(prob)) molkill(sim,mptr,mptr->list,-1);
3411 return CMDok; }
3412
3413
3414 /* cmdkillmolinsphere */
cmdkillmolinsphere(simptr sim,cmdptr cmd,char * line2)3415 enum CMDcode cmdkillmolinsphere(simptr sim,cmdptr cmd,char *line2) {
3416 int itct,i,*index;
3417 char nm[STRCHAR];
3418 moleculeptr mptr;
3419 static enum MolecState ms;
3420 static int inscan=0,s;
3421
3422 if(inscan) goto scanportion;
3423 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3424
3425 if(!sim->srfss) return CMDok;
3426 i=molstring2index1(sim,line2,&ms,&index);
3427 SCMDCHECK(i!=-1,"species is missing or cannot be read");
3428 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
3429 SCMDCHECK(i!=-3,"cannot read molecule state value");
3430 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
3431 SCMDCHECK(i!=-7,"error allocating memory");
3432 line2=strnword(line2,2);
3433 SCMDCHECK(line2,"missing surface name");
3434 itct=sscanf(line2,"%s",nm);
3435 SCMDCHECK(itct==1,"cannot read surface name");
3436 if(!strcmp(nm,"all")) s=-1;
3437 else {
3438 s=stringfind(sim->srfss->snames,sim->srfss->nsrf,nm);
3439 SCMDCHECK(s>=0,"surface not recognized"); }
3440
3441 if(i!=-4) {
3442 inscan=1;
3443 molscancmd(sim,i,index,ms,cmd,cmdkillmolinsphere);
3444 inscan=0; }
3445 return CMDok;
3446
3447 scanportion:
3448 mptr=(moleculeptr) line2;
3449 if(molinpanels(sim,mptr,s,PSsph))
3450 molkill(sim,mptr,mptr->list,-1);
3451 return CMDok; }
3452
3453
3454 /* cmdkillmolincmpt */
cmdkillmolincmpt(simptr sim,cmdptr cmd,char * line2)3455 enum CMDcode cmdkillmolincmpt(simptr sim,cmdptr cmd,char *line2) {
3456 int itct,i,c,*index;
3457 char cname[STRCHAR];
3458 moleculeptr mptr;
3459 enum MolecState ms;
3460 compartssptr cmptss;
3461 static compartptr cmpt;
3462 static int inscan=0;
3463
3464 if(inscan) goto scanportion;
3465 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3466
3467 cmptss=sim->cmptss;
3468 SCMDCHECK(cmptss,"no compartments defined");
3469 SCMDCHECK(sim->mols,"molecules are undefined");
3470 SCMDCHECK(line2,"missing argument");
3471 i=molstring2index1(sim,line2,&ms,&index);
3472 SCMDCHECK(i!=-1,"species is missing or cannot be read");
3473 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
3474 SCMDCHECK(i!=-3,"cannot read molecule state value");
3475 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
3476 SCMDCHECK(i!=-7,"error allocating memory");
3477 SCMDCHECK(line2=strnword(line2,2),"missing value argument");
3478 itct=sscanf(line2,"%s",cname);
3479 SCMDCHECK(itct==1,"cannot read compartment name");
3480 c=stringfind(cmptss->cnames,cmptss->ncmpt,cname);
3481 SCMDCHECK(c>=0,"compartment name not recognized");
3482 cmpt=cmptss->cmptlist[c];
3483
3484 if(i!=-4) {
3485 inscan=1;
3486 molscancmd(sim,i,index,ms,cmd,cmdkillmolincmpt);
3487 inscan=0; }
3488 return CMDok;
3489
3490 scanportion:
3491 mptr=(moleculeptr) line2;
3492 if(posincompart(sim,mptr->pos,cmpt,0))
3493 molkill(sim,mptr,mptr->list,-1);
3494 return CMDok; }
3495
3496
3497 /* cmdkillmoloutsidesystem */
cmdkillmoloutsidesystem(simptr sim,cmdptr cmd,char * line2)3498 enum CMDcode cmdkillmoloutsidesystem(simptr sim,cmdptr cmd,char *line2) {
3499 int i,*index;
3500 moleculeptr mptr;
3501 enum MolecState ms;
3502 static int inscan=0;
3503
3504 if(inscan) goto scanportion;
3505 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3506
3507 if(!sim->srfss) return CMDok;
3508 i=molstring2index1(sim,line2,&ms,&index);
3509 SCMDCHECK(i!=-1,"species is missing or cannot be read");
3510 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
3511 SCMDCHECK(i!=-3,"cannot read molecule state value");
3512 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
3513 SCMDCHECK(i!=-7,"error allocating memory");
3514
3515 if(i!=-4) {
3516 inscan=1;
3517 molscancmd(sim,i,index,ms,cmd,cmdkillmoloutsidesystem);
3518 inscan=0; }
3519 return CMDok;
3520
3521 scanportion:
3522 mptr=(moleculeptr) line2;
3523 if(!posinsystem(sim,mptr->pos)) molkill(sim,mptr,mptr->list,-1);
3524 return CMDok; }
3525
3526
3527 /* cmdfixmolcount */
cmdfixmolcount(simptr sim,cmdptr cmd,char * line2)3528 enum CMDcode cmdfixmolcount(simptr sim,cmdptr cmd,char *line2) {
3529 int itct,num,i,ll,m,ct,numl;
3530 char nm[STRCHAR];
3531 double pos1[DIMMAX],pos2[DIMMAX];
3532
3533 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3534
3535 SCMDCHECK(line2,"missing argument");
3536 SCMDCHECK(sim->mols,"molecules are undefined");
3537 itct=sscanf(line2,"%s %i",nm,&num);
3538 SCMDCHECK(itct==2,"read failure");
3539 SCMDCHECK(num>=0,"number cannot be negative");
3540 i=stringfind(sim->mols->spname,sim->mols->nspecies,nm);
3541 SCMDCHECK(i>=1,"name not recognized");
3542
3543 ll=sim->mols->listlookup[i][MSsoln];
3544 numl=sim->mols->nl[ll];
3545 ct=0;
3546 for(m=0;m<numl;m++)
3547 if(sim->mols->live[ll][m]->ident==i) ct++;
3548
3549 if(ct==num);
3550 else if(ct<num) {
3551 systemcorners(sim,pos1,pos2);
3552 SCMDCHECK(addmol(sim,num-ct,i,pos1,pos2,1)==0,"not enough available molecules"); }
3553 else {
3554 num=ct-num;
3555 for(;num>0;num--) {
3556 m=intrand(numl);
3557 while(sim->mols->live[ll][m]->ident!=i) m=(m==numl-1)?0:m+1;
3558 molkill(sim,sim->mols->live[ll][m],ll,m); }}
3559
3560 return CMDok; }
3561
3562
3563 /* cmdfixmolcountrange */
cmdfixmolcountrange(simptr sim,cmdptr cmd,char * line2)3564 enum CMDcode cmdfixmolcountrange(simptr sim,cmdptr cmd,char *line2) {
3565 int itct,lownum,highnum,i,ll,m,ct,numl;
3566 char nm[STRCHAR];
3567 double pos1[DIMMAX],pos2[DIMMAX];
3568
3569 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3570
3571 SCMDCHECK(line2,"missing argument");
3572 SCMDCHECK(sim->mols,"molecules are undefined");
3573 itct=strmathsscanf(line2,"%s %mi %mi",Varnames,Varvalues,Nvar,nm,&lownum,&highnum);
3574 SCMDCHECK(itct==3,"read failure");
3575 SCMDCHECK(lownum>=0 && highnum>=0 && highnum>=lownum,"molecule numbers are out of bounds");
3576 i=stringfind(sim->mols->spname,sim->mols->nspecies,nm);
3577 SCMDCHECK(i>=1,"species name not recognized");
3578
3579 ll=sim->mols->listlookup[i][MSsoln];
3580 numl=sim->mols->nl[ll];
3581 ct=0;
3582 for(m=0;m<numl;m++)
3583 if(sim->mols->live[ll][m]->ident==i) ct++;
3584
3585 if(ct>=lownum && ct<=highnum);
3586 else if(ct<lownum) {
3587 systemcorners(sim,pos1,pos2);
3588 SCMDCHECK(addmol(sim,lownum-ct,i,pos1,pos2,1)==0,"not enough available molecules"); }
3589 else {
3590 highnum=ct-highnum;
3591 for(;highnum>0;highnum--) {
3592 m=intrand(numl);
3593 while(sim->mols->live[ll][m]->ident!=i) m=(m==numl-1)?0:m+1;
3594 molkill(sim,sim->mols->live[ll][m],ll,m); }}
3595
3596 return CMDok; }
3597
3598
3599 /* cmdfixmolcountonsurf */
cmdfixmolcountonsurf(simptr sim,cmdptr cmd,char * line2)3600 enum CMDcode cmdfixmolcountonsurf(simptr sim,cmdptr cmd,char *line2) {
3601 int itct,num,i,ll,m,ct,numl,s,*index;
3602 char nm[STRCHAR];
3603 enum MolecState ms;
3604 surfaceptr sptr;
3605 moleculeptr mptr;
3606
3607 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3608 SCMDCHECK(line2,"missing argument");
3609 i=molstring2index1(sim,line2,&ms,&index);
3610 SCMDCHECK(i!=-1,"species is missing or cannot be read");
3611 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
3612 SCMDCHECK(i!=-3,"cannot read molecule state value");
3613 SCMDCHECK(i!=-4,"molecule name not recognized");
3614 SCMDCHECK(i!=-7,"error allocating memory");
3615 SCMDCHECK(i>0,"molecule name needs to be for a single species");
3616 SCMDCHECK(ms!=MSsoln && ms!=MSbsoln,"molecule state needs to be surface-bound");
3617 line2=strnword(line2,2);
3618 SCMDCHECK(line2,"fixmolcountonsurf format: species(state) number surface");
3619 itct=strmathsscanf(line2,"%mi %s",Varnames,Varvalues,Nvar,&num,nm);
3620 SCMDCHECK(itct==2,"read failure");
3621 SCMDCHECK(num>=0,"number cannot be negative");
3622 SCMDCHECK(sim->srfss,"no surfaces defined");
3623 s=stringfind(sim->srfss->snames,sim->srfss->nsrf,nm);
3624 SCMDCHECK(s>=0,"surface not recognized");
3625 sptr=sim->srfss->srflist[s];
3626
3627 ll=sim->mols->listlookup[i][ms];
3628 numl=sim->mols->nl[ll];
3629 ct=0;
3630 for(m=0;m<numl;m++) {
3631 mptr=sim->mols->live[ll][m];
3632 if(mptr->ident==i && mptr->mstate==ms && mptr->pnl->srf==sptr) ct++; }
3633
3634 if(ct==num);
3635 else if(ct<num) {
3636 SCMDCHECK(addsurfmol(sim,num-ct,i,ms,NULL,NULL,s,PSall,NULL)==0,"not enough available molecules"); }
3637 else {
3638 num=ct-num;
3639 for(;num>0;num--) {
3640 m=intrand(numl);
3641 mptr=sim->mols->live[ll][m];
3642 while(!(mptr->ident==i && mptr->mstate==ms && mptr->pnl->srf==sptr)) {
3643 m=(m==numl-1)?0:m+1;
3644 mptr=sim->mols->live[ll][m]; }
3645 molkill(sim,mptr,ll,m); }}
3646
3647 return CMDok; }
3648
3649
3650 /* cmdfixmolcountrangeonsurf */
cmdfixmolcountrangeonsurf(simptr sim,cmdptr cmd,char * line2)3651 enum CMDcode cmdfixmolcountrangeonsurf(simptr sim,cmdptr cmd,char *line2) {
3652 int itct,lownum,highnum,i,ll,m,ct,numl,s,*index;
3653 char nm[STRCHAR];
3654 enum MolecState ms;
3655 surfaceptr sptr;
3656 moleculeptr mptr;
3657
3658 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3659
3660 SCMDCHECK(line2,"missing argument");
3661 i=molstring2index1(sim,line2,&ms,&index);
3662 SCMDCHECK(i!=-1,"species is missing or cannot be read");
3663 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
3664 SCMDCHECK(i!=-3,"cannot read molecule state value");
3665 SCMDCHECK(i!=-4,"molecule name not recognized");
3666 SCMDCHECK(i!=-7,"error allocating memory");
3667 SCMDCHECK(i>0,"molecule name needs to be for a single species");
3668 SCMDCHECK(ms!=MSsoln && ms!=MSbsoln,"molecule state needs to be surface-bound");
3669 line2=strnword(line2,2);
3670 SCMDCHECK(line2,"fixmolcountrangeonsurf format: species(state) low_number high_number surface");
3671 itct=strmathsscanf(line2,"%mi %mi %s",Varnames,Varvalues,Nvar,&lownum,&highnum,nm);
3672 SCMDCHECK(itct==3,"read failure");
3673 SCMDCHECK(lownum>=0 && highnum>=0 && highnum>=lownum,"molecule numbers are out of bounds");
3674 SCMDCHECK(sim->srfss,"no surfaces defined");
3675 s=stringfind(sim->srfss->snames,sim->srfss->nsrf,nm);
3676 SCMDCHECK(s>=0,"surface not recognized");
3677 sptr=sim->srfss->srflist[s];
3678
3679 ll=sim->mols->listlookup[i][ms];
3680 numl=sim->mols->nl[ll];
3681 ct=0;
3682 for(m=0;m<numl;m++) {
3683 mptr=sim->mols->live[ll][m];
3684 if(mptr->ident==i && mptr->mstate==ms && mptr->pnl->srf==sptr) ct++; }
3685
3686 if(ct>=lownum && ct<=highnum);
3687 else if(ct<lownum) {
3688 SCMDCHECK(addsurfmol(sim,lownum-ct,i,ms,NULL,NULL,s,PSall,NULL)==0,"not enough available molecules"); }
3689 else {
3690 highnum=ct-highnum;
3691 for(;highnum>0;highnum--) {
3692 m=intrand(numl);
3693 mptr=sim->mols->live[ll][m];
3694 while(!(mptr->ident==i && mptr->mstate==ms && mptr->pnl->srf==sptr)) {
3695 m=(m==numl-1)?0:m+1;
3696 mptr=sim->mols->live[ll][m]; }
3697 molkill(sim,mptr,ll,m); }}
3698
3699 return CMDok; }
3700
3701
3702 /* cmdfixmolcountincmpt */
cmdfixmolcountincmpt(simptr sim,cmdptr cmd,char * line2)3703 enum CMDcode cmdfixmolcountincmpt(simptr sim,cmdptr cmd,char *line2) {
3704 int itct,num,i,ll,m,ct,numl,c;
3705 char nm[STRCHAR];
3706 moleculeptr mptr;
3707 compartptr cmpt;
3708
3709 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3710
3711 SCMDCHECK(line2,"missing argument");
3712 SCMDCHECK(sim->mols,"molecules are undefined");
3713 SCMDCHECK(sim->cmptss,"compartments are undefined");
3714 itct=strmathsscanf(line2,"%s %mi",Varnames,Varvalues,Nvar,nm,&num);
3715 SCMDCHECK(itct==2,"read failure");
3716 SCMDCHECK(num>=0,"number cannot be negative");
3717 i=stringfind(sim->mols->spname,sim->mols->nspecies,nm);
3718 SCMDCHECK(i>=1,"molecule name not recognized");
3719 line2=strnword(line2,3);
3720 SCMDCHECK(line2,"compartment name missing");
3721 itct=sscanf(line2,"%s",nm);
3722 c=stringfind(sim->cmptss->cnames,sim->cmptss->ncmpt,nm);
3723 SCMDCHECK(c>=0,"compartment not recognized");
3724 cmpt=sim->cmptss->cmptlist[c];
3725
3726 ll=sim->mols->listlookup[i][MSsoln];
3727 numl=sim->mols->nl[ll];
3728 ct=0;
3729 for(m=0;m<numl;m++) {
3730 mptr=sim->mols->live[ll][m];
3731 if(mptr->ident==i && mptr->mstate==MSsoln && posincompart(sim,mptr->pos,cmpt,0)) ct++; }
3732
3733 if(ct==num);
3734 else if(ct<num) {
3735 SCMDCHECK(addcompartmol(sim,num-ct,i,cmpt)==0,"not enough available molecules"); }
3736 else {
3737 num=ct-num;
3738 for(;num>0;num--) {
3739 m=intrand(numl);
3740 mptr=sim->mols->live[ll][m];
3741 while(!(mptr->ident==i && mptr->mstate==MSsoln && posincompart(sim,mptr->pos,cmpt,0))) {
3742 m=(m==numl-1)?0:m+1;
3743 mptr=sim->mols->live[ll][m]; }
3744 molkill(sim,mptr,ll,m); }}
3745
3746 return CMDok; }
3747
3748
3749 /* cmdfixmolcountrangeincmpt */
cmdfixmolcountrangeincmpt(simptr sim,cmdptr cmd,char * line2)3750 enum CMDcode cmdfixmolcountrangeincmpt(simptr sim,cmdptr cmd,char *line2) {
3751 int itct,i,ll,m,ct,numl,c,lownum,highnum;
3752 char nm[STRCHAR];
3753 moleculeptr mptr;
3754 compartptr cmpt;
3755
3756 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3757
3758 SCMDCHECK(line2,"missing argument");
3759 SCMDCHECK(sim->mols,"molecules are undefined");
3760 SCMDCHECK(sim->cmptss,"compartments are undefined");
3761 itct=strmathsscanf(line2,"%s %mi %mi",Varnames,Varvalues,Nvar,nm,&lownum,&highnum);
3762 SCMDCHECK(itct==3,"read failure");
3763 i=stringfind(sim->mols->spname,sim->mols->nspecies,nm);
3764 SCMDCHECK(i>=1,"molecule name not recognized");
3765 line2=strnword(line2,4);
3766 SCMDCHECK(line2,"compartment name missing");
3767 itct=sscanf(line2,"%s",nm);
3768 c=stringfind(sim->cmptss->cnames,sim->cmptss->ncmpt,nm);
3769 SCMDCHECK(c>=0,"compartment not recognized");
3770 cmpt=sim->cmptss->cmptlist[c];
3771
3772 ll=sim->mols->listlookup[i][MSsoln];
3773 numl=sim->mols->nl[ll];
3774 ct=0;
3775 for(m=0;m<numl;m++) {
3776 mptr=sim->mols->live[ll][m];
3777 if(mptr->ident==i && mptr->mstate==MSsoln && posincompart(sim,mptr->pos,cmpt,0)) ct++; }
3778
3779 if(ct>=lownum && ct<=highnum);
3780 else if(ct<lownum) {
3781 SCMDCHECK(addcompartmol(sim,lownum-ct,i,cmpt)==0,"not enough available molecules"); }
3782 else {
3783 highnum=ct-highnum;
3784 for(;highnum>0;highnum--) {
3785 m=intrand(numl);
3786 mptr=sim->mols->live[ll][m];
3787 while(!(mptr->ident==i && mptr->mstate==MSsoln && posincompart(sim,mptr->pos,cmpt,0))) {
3788 m=(m==numl-1)?0:m+1;
3789 mptr=sim->mols->live[ll][m]; }
3790 molkill(sim,mptr,ll,m); }}
3791
3792 return CMDok; }
3793
3794
3795 /* cmdequilmol */
cmdequilmol(simptr sim,cmdptr cmd,char * line2)3796 enum CMDcode cmdequilmol(simptr sim,cmdptr cmd,char *line2) {
3797 int itct,*index;
3798 moleculeptr mptr;
3799 static enum MolecState ms1,ms2;
3800 static double prob;
3801 static int xyzvar;
3802 static char probstr[STRCHAR];
3803 static int inscan=0,i1,i2;
3804
3805 if(inscan) goto scanportion;
3806 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3807
3808 i1=molstring2index1(sim,line2,&ms1,&index);
3809 SCMDCHECK(i1!=-1,"species is missing or cannot be read");
3810 SCMDCHECK(i1!=-2,"mismatched or improper parentheses around molecule state");
3811 SCMDCHECK(i1!=-3,"cannot read molecule state value");
3812 SCMDCHECK(i1!=-4,"molecule name not recognized");
3813 SCMDCHECK(i1!=-7,"error allocating memory");
3814 SCMDCHECK(i1>0,"molecule name has to be for a single species");
3815 SCMDCHECK(ms1!=MSall,"molecule state cannot be 'all'");
3816 line2=strnword(line2,2);
3817 i2=molstring2index1(sim,line2,&ms2,&index);
3818 SCMDCHECK(i2!=-1,"species is missing or cannot be read");
3819 SCMDCHECK(i2!=-2,"mismatched or improper parentheses around molecule state");
3820 SCMDCHECK(i2!=-3,"cannot read molecule state value");
3821 SCMDCHECK(i2!=-4,"molecule name not recognized");
3822 SCMDCHECK(i2!=-7,"error allocating memory");
3823 SCMDCHECK(i2>0,"molecule name has to be for a single species");
3824 SCMDCHECK(ms2!=MSall,"molecule state cannot be 'all'");
3825 SCMDCHECK((ms1==MSsoln && ms2==MSsoln) || (ms1!=MSsoln && ms2!=MSsoln),"cannot equilibrate between solution and surface-bound");
3826 line2=strnword(line2,2);
3827 SCMDCHECK(line2,"missing probability argument");
3828 if(strhasname(line2,"x") || strhasname(line2,"y") || strhasname(line2,"z")) {
3829 xyzvar=1;
3830 strcpy(probstr,line2); }
3831 else {
3832 xyzvar=0;
3833 itct=strmathsscanf(line2,"%mlg",Varnames,Varvalues,Nvar,&prob);
3834 SCMDCHECK(itct==1,"failed to read probability");
3835 SCMDCHECK(prob>=0 && prob<=1,"probability is out of bounds"); }
3836
3837 inscan=1;
3838 molscancmd(sim,-1,index,MSall,cmd,cmdequilmol);
3839 inscan=0;
3840 return CMDok;
3841
3842 scanportion:
3843 mptr=(moleculeptr) line2;
3844 if((mptr->ident==i1 && mptr->mstate==ms1) || (mptr->ident==i2 && mptr->mstate==ms2)) {
3845 if(xyzvar) {
3846 simsetvariable(sim,"x",mptr->pos[0]);
3847 if(sim->dim>1) simsetvariable(sim,"y",mptr->pos[1]);
3848 if(sim->dim>2) simsetvariable(sim,"z",mptr->pos[2]);
3849 strmathsscanf(probstr,"%mlg",Varnames,Varvalues,Nvar,&prob); }
3850 if(coinrandD(prob))
3851 molchangeident(sim,mptr,-1,-1,i2,ms2,mptr->pnl);
3852 else
3853 molchangeident(sim,mptr,-1,-1,i1,ms1,mptr->pnl); }
3854 return CMDok; }
3855
3856
3857 /* cmdreplacemol */
cmdreplacemol(simptr sim,cmdptr cmd,char * line2)3858 enum CMDcode cmdreplacemol(simptr sim,cmdptr cmd,char *line2) {
3859 int itct,i1,*index1,*index2;
3860 enum MolecState ms1;
3861 moleculeptr mptr;
3862 static enum MolecState ms2;
3863 static double prob;
3864 static char probstr[STRCHAR];
3865 static int xyzvar;
3866 static int inscan=0,i2;
3867
3868 if(inscan) goto scanportion;
3869 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3870
3871 i1=molstring2index1(sim,line2,&ms1,&index1);
3872 SCMDCHECK(i1!=-1,"species is missing or cannot be read");
3873 SCMDCHECK(i1!=-2,"mismatched or improper parentheses around molecule state");
3874 SCMDCHECK(i1!=-3,"cannot read molecule state value");
3875 SCMDCHECK(i1!=-4,"molecule name not recognized");
3876 SCMDCHECK(i1!=-7,"error allocating memory");
3877 SCMDCHECK(ms1!=MSall,"molecule state cannot be 'all'");
3878 line2=strnword(line2,2);
3879 SCMDCHECK(line2,"missing second species name");
3880 i2=molstring2index1(sim,line2,&ms2,&index2);
3881 SCMDCHECK(i2!=-1,"species is missing or cannot be read");
3882 SCMDCHECK(i2!=-2,"mismatched or improper parentheses around molecule state");
3883 SCMDCHECK(i2!=-3,"cannot read molecule state value");
3884 SCMDCHECK(i2!=-4,"molecule name not recognized");
3885 SCMDCHECK(i2!=-7,"error allocating memory");
3886 SCMDCHECK(i2>0,"molecule name has to be for a single species");
3887 SCMDCHECK(ms2!=MSall,"molecule state cannot be 'all'");
3888 SCMDCHECK((ms1==MSsoln && ms2==MSsoln) || (ms1!=MSsoln && ms2!=MSsoln),"cannot equilibrate between solution and surface-bound");
3889 line2=strnword(line2,2);
3890 SCMDCHECK(line2,"missing probability information");
3891 itct=sscanf(line2,"%s",probstr);
3892 SCMDCHECK(itct==1,"missing probability information");
3893 if(strhasname(probstr,"x") || strhasname(probstr,"y") || strhasname(probstr,"z"))
3894 xyzvar=1;
3895 else {
3896 xyzvar=0;
3897 itct=strmathsscanf(line2,"%mlg",Varnames,Varvalues,Nvar,&prob);
3898 SCMDCHECK(itct==1,"cannot read fraction");
3899 SCMDCHECK(prob>=0 && prob<=1,"fraction out of bounds"); }
3900
3901 inscan=1;
3902 molscancmd(sim,i1,index1,ms1,cmd,cmdreplacemol);
3903 inscan=0;
3904 return CMDok;
3905
3906 scanportion:
3907 mptr=(moleculeptr) line2;
3908 if(xyzvar) {
3909 simsetvariable(sim,"x",mptr->pos[0]);
3910 if(sim->dim>1) simsetvariable(sim,"y",mptr->pos[1]);
3911 if(sim->dim>2) simsetvariable(sim,"z",mptr->pos[2]);
3912 strmathsscanf(probstr,"%mlg",Varnames,Varvalues,Nvar,&prob); }
3913 if(coinrandD(prob))
3914 molchangeident(sim,mptr,-1,-1,i2,ms2,mptr->pnl);
3915 return CMDok; }
3916
3917
3918 /* cmdreplacexyzmol */
cmdreplacexyzmol(simptr sim,cmdptr cmd,char * line2)3919 enum CMDcode cmdreplacexyzmol(simptr sim,cmdptr cmd,char *line2) {
3920 int itct,i,m,d,ll,*index;
3921 moleculeptr *mlist;
3922 boxptr bptr;
3923 double pos[DIMMAX];
3924 enum MolecState ms;
3925
3926 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3927
3928 i=molstring2index1(sim,line2,&ms,&index);
3929 SCMDCHECK(i!=-1,"species is missing or cannot be read");
3930 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
3931 SCMDCHECK(i!=-3,"cannot read molecule state value");
3932 SCMDCHECK(i!=-4,"molecule name not recognized");
3933 SCMDCHECK(i!=-7,"error allocating memory");
3934 SCMDCHECK(i>0,"molecule name has to be for a single species");
3935 SCMDCHECK(ms!=MSall,"molecule state cannot be 'all'");
3936 line2=strnword(line2,2);
3937 SCMDCHECK(line2,"missing position information");
3938 if(sim->dim==1) itct=strmathsscanf(line2,"%mlg",Varnames,Varvalues,Nvar,&pos[0]);
3939 else if(sim->dim==2) itct=strmathsscanf(line2,"%mlg %mlg",Varnames,Varvalues,Nvar,&pos[0],&pos[1]);
3940 else itct=strmathsscanf(line2,"%mlg %mlg %mlg",Varnames,Varvalues,Nvar,&pos[0],&pos[1],&pos[2]);
3941 SCMDCHECK(itct==sim->dim,"insufficient dimensions entered");
3942 bptr=pos2box(sim,pos);
3943 ll=sim->mols->listlookup[i][ms];
3944 mlist=bptr->mol[ll];
3945 for(m=0;m<bptr->nmol[ll];m++) {
3946 for(d=0;d<sim->dim;d++)
3947 if(mlist[m]->pos[d]!=pos[d]) d=sim->dim+1;
3948 if(d==sim->dim) {
3949 molchangeident(sim,mlist[m],ll,-1,i,ms,mlist[m]->pnl);
3950 m=bptr->nmol[ll]+1; }}
3951 return CMDok; }
3952
3953
3954 /* cmdreplacevolmol */
cmdreplacevolmol(simptr sim,cmdptr cmd,char * line2)3955 enum CMDcode cmdreplacevolmol(simptr sim,cmdptr cmd,char *line2) {
3956 int m,itct,dim,d,b,b1,b2,i1,i2,ll,*index;
3957 double *pos,poslo[DIMMAX],poshi[DIMMAX],frac;
3958 boxptr bptr1,bptr2,bptr;
3959 boxssptr boxs;
3960 moleculeptr *mlist,mptr;
3961 enum MolecState ms1,ms2;
3962 char probstr[STRCHAR];
3963 int xyzvar;
3964
3965 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
3966
3967 i1=molstring2index1(sim,line2,&ms1,&index);
3968 SCMDCHECK(i1!=-1,"species is missing or cannot be read");
3969 SCMDCHECK(i1!=-2,"mismatched or improper parentheses around molecule state");
3970 SCMDCHECK(i1!=-3,"cannot read molecule state value");
3971 SCMDCHECK(i1!=-4,"molecule name not recognized");
3972 SCMDCHECK(i1!=-7,"error allocating memory");
3973 SCMDCHECK(i1>0,"molecule name has to be for a single species");
3974 SCMDCHECK(ms1!=MSall,"molecule state cannot be 'all'");
3975 line2=strnword(line2,2);
3976 i2=molstring2index1(sim,line2,&ms2,&index);
3977 SCMDCHECK(i2!=-1,"species is missing or cannot be read");
3978 SCMDCHECK(i2!=-2,"mismatched or improper parentheses around molecule state");
3979 SCMDCHECK(i2!=-3,"cannot read molecule state value");
3980 SCMDCHECK(i2!=-4,"molecule name not recognized");
3981 SCMDCHECK(i2!=-7,"error allocating memory");
3982 SCMDCHECK(i2>0,"molecule name has to be for a single species");
3983 SCMDCHECK(ms2!=MSall,"molecule state cannot be 'all'");
3984 SCMDCHECK((ms1==MSsoln && ms2==MSsoln) || (ms1!=MSsoln && ms2!=MSsoln),"cannot equilibrate between solution and surface-bound");
3985 line2=strnword(line2,2);
3986 SCMDCHECK(line2,"missing fraction information");
3987 itct=sscanf(line2,"%s",probstr);
3988 SCMDCHECK(itct==1,"cannot read fraction");
3989 if(strhasname(probstr,"x") || strhasname(probstr,"y") || strhasname(probstr,"z"))
3990 xyzvar=1;
3991 else {
3992 xyzvar=0;
3993 itct=strmathsscanf(line2,"%mlg",Varnames,Varvalues,Nvar,&frac);
3994 SCMDCHECK(itct==1,"cannot read fraction");
3995 SCMDCHECK(frac>=0 && frac<=1,"fraction out of bounds"); }
3996 line2=strnword(line2,2);
3997 dim=sim->dim;
3998 boxs=sim->boxs;
3999 for(d=0;d<dim;d++) {
4000 SCMDCHECK(line2,"missing argument");
4001 itct=strmathsscanf(line2,"%mlg %mlg",Varnames,Varvalues,Nvar,&poslo[d],&poshi[d]);
4002 SCMDCHECK(itct==2,"read failure");
4003 line2=strnword(line2,3); }
4004
4005 bptr1=pos2box(sim,poslo);
4006 bptr2=pos2box(sim,poshi);
4007 b1=indx2addZV(bptr1->indx,boxs->side,dim);
4008 b2=indx2addZV(bptr2->indx,boxs->side,dim);
4009 ll=sim->mols->listlookup[i1][ms1];
4010 for(b=b1;b<=b2;b=nextaddZV(b,bptr1->indx,bptr2->indx,boxs->side,dim)) {
4011 bptr=boxs->blist[b];
4012 mlist=bptr->mol[ll];
4013 for(m=0;m<bptr->nmol[ll];m++) {
4014 mptr=mlist[m];
4015 if(mptr->ident==i1 && mptr->mstate==ms1) {
4016 pos=mptr->pos;
4017 for(d=0;d<dim;d++) if(pos[d]<poslo[d] || pos[d]>poshi[d]) d=dim+1;
4018 if(d==dim) {
4019 if(xyzvar) {
4020 simsetvariable(sim,"x",mptr->pos[0]);
4021 if(sim->dim>1) simsetvariable(sim,"y",mptr->pos[1]);
4022 if(sim->dim>2) simsetvariable(sim,"z",mptr->pos[2]);
4023 strmathsscanf(probstr,"%mlg",Varnames,Varvalues,Nvar,&frac); }
4024 if(coinrandD(frac)) {
4025 molchangeident(sim,mptr,ll,-1,i2,ms2,mptr->pnl); }}}}}
4026 return CMDok; }
4027
4028
4029 /* cmdreplacecmptmol */
cmdreplacecmptmol(simptr sim,cmdptr cmd,char * line2)4030 enum CMDcode cmdreplacecmptmol(simptr sim,cmdptr cmd,char *line2) {
4031 int itct,i1,c,*index1,*index2;
4032 char nm[STRCHAR];
4033 enum MolecState ms1;
4034 moleculeptr mptr;
4035 static enum MolecState ms2;
4036 static compartptr cmpt;
4037 static double frac;
4038 static char probstr[STRCHAR];
4039 static int xyzvar;
4040 static int inscan=0,i2;
4041
4042 if(inscan) goto scanportion;
4043 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
4044
4045 i1=molstring2index1(sim,line2,&ms1,&index1);
4046 SCMDCHECK(i1!=-1,"species is missing or cannot be read");
4047 SCMDCHECK(i1!=-2,"mismatched or improper parentheses around molecule state");
4048 SCMDCHECK(i1!=-3,"cannot read molecule state value");
4049 SCMDCHECK(i1!=-4,"molecule name not recognized");
4050 SCMDCHECK(i1!=-7,"error allocating memory");
4051 SCMDCHECK(ms1!=MSall,"molecule state cannot be 'all'");
4052 line2=strnword(line2,2);
4053 SCMDCHECK(line2,"missing second species name");
4054 i2=molstring2index1(sim,line2,&ms2,&index2);
4055 SCMDCHECK(i2!=-1,"species is missing or cannot be read");
4056 SCMDCHECK(i2!=-2,"mismatched or improper parentheses around molecule state");
4057 SCMDCHECK(i2!=-3,"cannot read molecule state value");
4058 SCMDCHECK(i2!=-4,"molecule name not recognized");
4059 SCMDCHECK(i2!=-7,"error allocating memory");
4060 SCMDCHECK(i2>0,"molecule name has to be for a single species");
4061 SCMDCHECK(ms2!=MSall,"molecule state cannot be 'all'");
4062 SCMDCHECK((ms1==MSsoln && ms2==MSsoln) || (ms1!=MSsoln && ms2!=MSsoln),"cannot equilibrate between solution and surface-bound");
4063 line2=strnword(line2,2);
4064 SCMDCHECK(line2,"missing fraction information");
4065 itct=sscanf(line2,"%s",probstr);
4066 SCMDCHECK(itct==1,"missing fraction information");
4067 if(strhasname(probstr,"x") || strhasname(probstr,"y") || strhasname(probstr,"z"))
4068 xyzvar=1;
4069 else {
4070 xyzvar=0;
4071 itct=strmathsscanf(line2,"%mlg",Varnames,Varvalues,Nvar,&frac);
4072 SCMDCHECK(itct==1,"cannot read fraction");
4073 SCMDCHECK(frac>=0 && frac<=1,"fraction out of bounds"); }
4074 line2=strnword(line2,2);
4075 SCMDCHECK(line2,"compartment name missing");
4076 itct=sscanf(line2,"%s",nm);
4077 c=stringfind(sim->cmptss->cnames,sim->cmptss->ncmpt,nm);
4078 SCMDCHECK(c>=0,"compartment not recognized");
4079 cmpt=sim->cmptss->cmptlist[c];
4080
4081 inscan=1;
4082 molscancmd(sim,i1,index1,ms1,cmd,cmdreplacecmptmol);
4083 inscan=0;
4084 return CMDok;
4085
4086 scanportion:
4087 mptr=(moleculeptr) line2;
4088 if(posincompart(sim,mptr->pos,cmpt,0)) {
4089 if(xyzvar) {
4090 simsetvariable(sim,"x",mptr->pos[0]);
4091 if(sim->dim>1) simsetvariable(sim,"y",mptr->pos[1]);
4092 if(sim->dim>2) simsetvariable(sim,"z",mptr->pos[2]);
4093 strmathsscanf(probstr,"%mlg",Varnames,Varvalues,Nvar,&frac); }
4094 if(coinrandD(frac))
4095 molchangeident(sim,mptr,-1,-1,i2,ms2,mptr->pnl); }
4096 return CMDok; }
4097
4098
4099 /* cmdmodulatemol */
cmdmodulatemol(simptr sim,cmdptr cmd,char * line2)4100 enum CMDcode cmdmodulatemol(simptr sim,cmdptr cmd,char *line2) {
4101 int itct,*index;
4102 moleculeptr mptr;
4103 double freq,shift;
4104 static double prob;
4105 static enum MolecState ms1,ms2;
4106 static int inscan=0,i1,i2;
4107
4108 if(inscan) goto scanportion;
4109 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
4110
4111 i1=molstring2index1(sim,line2,&ms1,&index);
4112 SCMDCHECK(i1!=-1,"species is missing or cannot be read");
4113 SCMDCHECK(i1!=-2,"mismatched or improper parentheses around molecule state");
4114 SCMDCHECK(i1!=-3,"cannot read molecule state value");
4115 SCMDCHECK(i1!=-4,"molecule name not recognized");
4116 SCMDCHECK(i1!=-7,"error allocating memory");
4117 SCMDCHECK(i1>0,"molecule name has to be for a single species");
4118 SCMDCHECK(ms1!=MSall,"molecule state cannot be 'all'");
4119 line2=strnword(line2,2);
4120 i2=molstring2index1(sim,line2,&ms2,&index);
4121 SCMDCHECK(i2!=-1,"species is missing or cannot be read");
4122 SCMDCHECK(i2!=-2,"mismatched or improper parentheses around molecule state");
4123 SCMDCHECK(i2!=-3,"cannot read molecule state value");
4124 SCMDCHECK(i2!=-4,"molecule name not recognized");
4125 SCMDCHECK(i2!=-7,"error allocating memory");
4126 SCMDCHECK(i2>0,"molecule name has to be for a single species");
4127 SCMDCHECK(ms2!=MSall,"molecule state cannot be 'all'");
4128 SCMDCHECK((ms1==MSsoln && ms2==MSsoln) || (ms1!=MSsoln && ms2!=MSsoln),"cannot equilibrate between solution and surface-bound");
4129 line2=strnword(line2,2);
4130 SCMDCHECK(line2,"missing frequency and shift");
4131 itct=strmathsscanf(line2,"%mlg %mlg",Varnames,Varvalues,Nvar,&freq,&shift);
4132 SCMDCHECK(itct==2,"failure reading frequency or shift");
4133
4134 inscan=1;
4135 molscancmd(sim,-1,index,MSall,cmd,cmdmodulatemol);
4136 inscan=0;
4137
4138 return CMDok;
4139
4140 scanportion:
4141 mptr=(moleculeptr) line2;
4142 if((mptr->ident==i1 && mptr->mstate==ms1) || (mptr->ident==i2 && mptr->mstate==ms2)) {
4143 if(coinrandD(prob))
4144 molchangeident(sim,mptr,-1,-1,i2,ms2,mptr->pnl);
4145 else
4146 molchangeident(sim,mptr,-1,-1,i1,ms1,mptr->pnl); }
4147 return CMDok; }
4148
4149
4150 /* cmdreact1 */
cmdreact1(simptr sim,cmdptr cmd,char * line2)4151 enum CMDcode cmdreact1(simptr sim,cmdptr cmd,char *line2) {
4152 int itct,i,r,*index;
4153 char rnm[STRCHAR];
4154 moleculeptr mptr;
4155 enum MolecState ms;
4156 static rxnptr rxn;
4157 static int inscan=0;
4158
4159 if(inscan) goto scanportion;
4160 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
4161
4162 i=molstring2index1(sim,line2,&ms,&index);
4163 SCMDCHECK(i!=-1,"species is missing or cannot be read");
4164 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
4165 SCMDCHECK(i!=-3,"cannot read molecule state value");
4166 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
4167 SCMDCHECK(i!=-7,"error allocating memory");
4168 SCMDCHECK(ms!=MSall,"molecule state cannot be 'all'");
4169 line2=strnword(line2,2);
4170 SCMDCHECK(line2,"reaction name is missing");
4171 itct=sscanf(line2,"%s",rnm);
4172 SCMDCHECK(itct==1,"cannot read reaction name");
4173 SCMDCHECK(sim->rxnss[1],"no first order reactions defined");
4174 r=stringfind(sim->rxnss[1]->rname,sim->rxnss[1]->totrxn,rnm);
4175 SCMDCHECK(r>=0,"reaction not recognized");
4176 rxn=sim->rxnss[1]->rxn[r];
4177
4178 if(i!=-4) {
4179 inscan=1;
4180 molscancmd(sim,i,index,ms,cmd,cmdreact1);
4181 inscan=0; }
4182 return CMDok;
4183
4184 scanportion:
4185 mptr=(moleculeptr) line2;
4186 doreact(sim,rxn,mptr,NULL,-1,-1,-1,-1,NULL,NULL);
4187 return CMDok; }
4188
4189
4190 /* cmdsetrateint */
cmdsetrateint(simptr sim,cmdptr cmd,char * line2)4191 enum CMDcode cmdsetrateint(simptr sim,cmdptr cmd,char *line2) {
4192 int itct,r,order;
4193 char rnm[STRCHAR];
4194 double rateint;
4195
4196 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
4197
4198 SCMDCHECK(line2,"missing argument");
4199 itct=strmathsscanf(line2,"%s %mlg",Varnames,Varvalues,Nvar,rnm,&rateint);
4200 SCMDCHECK(itct==2,"read failure");
4201 r=-1;
4202 if(sim->rxnss[0]) r=stringfind(sim->rxnss[0]->rname,sim->rxnss[0]->totrxn,rnm);
4203 if(r>=0) order=0;
4204 else {
4205 if(sim->rxnss[1]) r=stringfind(sim->rxnss[1]->rname,sim->rxnss[1]->totrxn,rnm);
4206 if(r>=0) order=1;
4207 else {
4208 if(sim->rxnss[2]) r=stringfind(sim->rxnss[2]->rname,sim->rxnss[2]->totrxn,rnm);
4209 if(r>=0) order=2;
4210 else SCMDCHECK(0,"reaction name not recognized"); }}
4211 SCMDCHECK(rateint>=0,"internal rate cannot be negative");
4212 if(order<2) RxnSetValue(sim,"prob",sim->rxnss[order]->rxn[r],rateint);
4213 else RxnSetValue(sim,"bindrad",sim->rxnss[order]->rxn[r],rateint);
4214 return CMDok; }
4215
4216
4217 /* cmdshufflemollist */
cmdshufflemollist(simptr sim,cmdptr cmd,char * line2)4218 enum CMDcode cmdshufflemollist(simptr sim,cmdptr cmd,char *line2) {
4219 int itct,ll,lllo,llhi;
4220 char nm[STRCHAR];
4221
4222 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
4223
4224 SCMDCHECK(line2,"missing argument");
4225 itct=sscanf(line2,"%s",nm);
4226 SCMDCHECK(itct==1,"read failure");
4227 SCMDCHECK(sim->mols && sim->mols->nlist>0,"no molecule lists");
4228 if(!strcmp(nm,"all")) ll=-1;
4229 else {
4230 ll=stringfind(sim->mols->listname,sim->mols->nlist,nm);
4231 SCMDCHECK(ll>=0,"list name not recognized"); }
4232 lllo=(ll==-1)?0:ll;
4233 llhi=(ll==-1)?sim->mols->nlist:ll+1;
4234 for(ll=lllo;ll<llhi;ll++)
4235 randshuffletableV((void**) sim->mols->live[ll],sim->mols->nl[ll]);
4236 return CMDok; }
4237
4238
4239 /* cmdshufflereactions */
cmdshufflereactions(simptr sim,cmdptr cmd,char * line2)4240 enum CMDcode cmdshufflereactions(simptr sim,cmdptr cmd,char *line2) {
4241 int itct,i,j,i1,i2,isym,j1,j2,*index1,*index2;
4242 char nm1[STRCHAR],nm2[STRCHAR];
4243 enum MolecState ms1,ms2;
4244 rxnssptr rxnss;
4245
4246 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
4247
4248 SCMDCHECK(line2,"missing argument");
4249 itct=sscanf(line2,"%s %s",nm1,nm2);
4250 SCMDCHECK(itct==2,"missing argument");
4251 i1=molstring2index1(sim,nm1,&ms1,&index1);
4252 SCMDCHECK(i1>=0 || i1==-5,"first species name not recognized");
4253 i2=molstring2index1(sim,nm2,&ms2,&index2);
4254 SCMDCHECK(i2>=0 || i2==-5,"second species name not recognized");
4255 rxnss=sim->rxnss[2];
4256 if(!rxnss) return CMDok;
4257
4258 for(j1=0;j1<index1[PDnresults];j1++)
4259 for(j2=0;j2<index2[PDnresults];j2++) {
4260 i1=index1[PDMAX+j1];
4261 i2=index2[PDMAX+j2];
4262 i=i1*rxnss->maxspecies+i2;
4263 if(rxnss->nrxn[i]) {
4264 randshuffletableI(rxnss->table[i],rxnss->nrxn[i]);
4265 isym=i2*rxnss->maxspecies+i1;
4266 for(j=0;j<rxnss->nrxn[i];j++)
4267 rxnss->table[isym][j]=rxnss->table[i][j]; }}
4268
4269 return CMDok; }
4270
4271
4272 /* cmdsettimestep */
cmdsettimestep(simptr sim,cmdptr cmd,char * line2)4273 enum CMDcode cmdsettimestep(simptr sim,cmdptr cmd,char *line2) {
4274 int itct,er;
4275 double dt;
4276
4277 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
4278
4279 SCMDCHECK(line2,"missing argument");
4280 itct=strmathsscanf(line2,"%mlg",Varnames,Varvalues,Nvar,&dt);
4281 SCMDCHECK(itct==1,"read failure");
4282 SCMDCHECK(dt>0,"time step must be >0");
4283
4284 er=simsettime(sim,dt,3);
4285 SCMDCHECK(er==0,"error while setting the simulation time step");
4286 return CMDok; }
4287
4288
4289 /* cmdporttransport */
cmdporttransport(simptr sim,cmdptr cmd,char * line2)4290 enum CMDcode cmdporttransport(simptr sim,cmdptr cmd,char *line2) {
4291 int itct,prt1,prt2;
4292 char nm1[STRCHAR],nm2[STRCHAR];
4293 portptr port1,port2;
4294
4295 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
4296
4297 SCMDCHECK(line2,"missing argument");
4298 SCMDCHECK(sim->portss,"no port superstructure is defined");
4299 itct=sscanf(line2,"%s %s",nm1,nm2);
4300 SCMDCHECK(itct==2,"porttransport format: port1 port2");
4301 prt1=stringfind(sim->portss->portnames,sim->portss->nport,nm1);
4302 SCMDCHECK(prt1>=0,"name of port1 is not recognized");
4303 prt2=stringfind(sim->portss->portnames,sim->portss->nport,nm2);
4304 SCMDCHECK(prt2>=0,"name of port2 is not recognized");
4305 port1=sim->portss->portlist[prt1];
4306 port2=sim->portss->portlist[prt2];
4307 porttransport(sim,port1,sim,port2);
4308 return CMDok; }
4309
4310
4311 /* cmdexcludebox */
cmdexcludebox(simptr sim,cmdptr cmd,char * line2)4312 enum CMDcode cmdexcludebox(simptr sim,cmdptr cmd,char *line2) {
4313 int m,itct,dim,d,b,b1,b2;
4314 double *pos,poslo[DIMMAX],poshi[DIMMAX];
4315 boxptr bptr1,bptr2,bptr;
4316 boxssptr boxs;
4317 moleculeptr *mlist;
4318
4319 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
4320
4321 dim=sim->dim;
4322 boxs=sim->boxs;
4323 for(d=0;d<dim;d++) {
4324 SCMDCHECK(line2,"missing argument");
4325 itct=strmathsscanf(line2,"%mlg %mlg",Varnames,Varvalues,Nvar,&poslo[d],&poshi[d]);
4326 SCMDCHECK(itct==2,"read failure");
4327 line2=strnword(line2,3); }
4328
4329 bptr1=pos2box(sim,poslo);
4330 bptr2=pos2box(sim,poshi);
4331 b1=indx2addZV(bptr1->indx,boxs->side,dim);
4332 b2=indx2addZV(bptr2->indx,boxs->side,dim);
4333 for(b=b1;b<=b2;b=nextaddZV(b,bptr1->indx,bptr2->indx,boxs->side,dim)) {
4334 bptr=boxs->blist[b];
4335 mlist=bptr->mol[0];
4336 for(m=0;m<bptr->nmol[0];m++) {
4337 pos=mlist[m]->pos;
4338 for(d=0;d<dim;d++) if(pos[d]<poslo[d] || pos[d]>poshi[d]) d=dim+1;
4339 if(d==dim) {
4340 pos=mlist[m]->posx;
4341 for(d=0;d<dim;d++) if(pos[d]<poslo[d] || pos[d]>poshi[d]) d=dim+1;
4342 if(d>dim) copyVD(mlist[m]->posx,mlist[m]->pos,dim); }}}
4343 sim->mols->touch++;
4344 return CMDok; }
4345
4346
4347 /* cmdexcludesphere */
cmdexcludesphere(simptr sim,cmdptr cmd,char * line2)4348 enum CMDcode cmdexcludesphere(simptr sim,cmdptr cmd,char *line2) {
4349 int m,itct,dim,d,b,b1,b2;
4350 double *pos,poslo[DIMMAX],poshi[DIMMAX],poscent[DIMMAX],rad,dist;
4351 boxptr bptr1,bptr2,bptr;
4352 boxssptr boxs;
4353 moleculeptr *mlist;
4354
4355 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
4356
4357 dim=sim->dim;
4358 boxs=sim->boxs;
4359 for(d=0;d<dim;d++) {
4360 SCMDCHECK(line2,"missing center argument");
4361 itct=strmathsscanf(line2,"%mlg",Varnames,Varvalues,Nvar,&poscent[d]);
4362 SCMDCHECK(itct==1,"failure reading center");
4363 line2=strnword(line2,2); }
4364 SCMDCHECK(line2,"missing radius");
4365 itct=strmathsscanf(line2,"%mlg",Varnames,Varvalues,Nvar,&rad);
4366 SCMDCHECK(itct==1,"failure reading radius");
4367
4368 dist=rad*sqrt((double)dim);
4369 for(d=0;d<dim;d++) {
4370 poslo[d]=poscent[d]-dist;
4371 poshi[d]=poscent[d]+dist; }
4372 rad*=rad;
4373 bptr1=pos2box(sim,poslo);
4374 bptr2=pos2box(sim,poshi);
4375 b1=indx2addZV(bptr1->indx,boxs->side,dim);
4376 b2=indx2addZV(bptr2->indx,boxs->side,dim);
4377 for(b=b1;b<=b2;b=nextaddZV(b,bptr1->indx,bptr2->indx,boxs->side,dim)) {
4378 bptr=boxs->blist[b];
4379 mlist=bptr->mol[0];
4380 for(m=0;m<bptr->nmol[0];m++) {
4381 pos=mlist[m]->pos;
4382 for(dist=0,d=0;d<dim;d++) if((dist+=(pos[d]-poscent[d])*(pos[d]-poscent[d]))>rad) d=dim+1;
4383 if(d==dim) {
4384 pos=mlist[m]->posx;
4385 for(dist=0,d=0;d<dim;d++) if((dist+=(pos[d]-poscent[d])*(pos[d]-poscent[d]))>rad) d=dim+1;
4386 if(d>dim) copyVD(mlist[m]->posx,mlist[m]->pos,dim); }}}
4387 sim->mols->touch++;
4388 return CMDok; }
4389
4390
4391 /* cmdincludeecoli */
cmdincludeecoli(simptr sim,cmdptr cmd,char * line2)4392 enum CMDcode cmdincludeecoli(simptr sim,cmdptr cmd,char *line2) {
4393 moleculeptr mptr;
4394 wallptr *wlist;
4395 static double rad,length,pos[DIMMAX];
4396 static int inscan=0;
4397
4398 if(inscan) goto scanportion;
4399 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
4400
4401 SCMDCHECK(sim->dim==3,"system is not 3 dimensional");
4402 wlist=sim->wlist;
4403 rad=0.5*(wlist[3]->pos-wlist[2]->pos);
4404 length=wlist[1]->pos-wlist[0]->pos;
4405 pos[0]=wlist[0]->pos;
4406 pos[1]=0.5*(wlist[2]->pos+wlist[3]->pos);
4407 pos[2]=0.5*(wlist[4]->pos+wlist[5]->pos);
4408
4409 inscan=1;
4410 molscancmd(sim,-1,NULL,MSsoln,cmd,cmdincludeecoli);
4411 inscan=0;
4412 sim->mols->touch++;
4413 return CMDok;
4414
4415 scanportion:
4416 mptr=(moleculeptr) line2;
4417 if(!insideecoli(mptr->pos,pos,rad,length)) {
4418 if(insideecoli(mptr->posx,pos,rad,length))
4419 copyVD(mptr->posx,mptr->pos,3);
4420 else putinecoli(mptr->pos,pos,rad,length); }
4421 return CMDok; }
4422
4423
4424 /* cmdsetreactionratemolcount */
cmdsetreactionratemolcount(simptr sim,cmdptr cmd,char * line2)4425 enum CMDcode cmdsetreactionratemolcount(simptr sim,cmdptr cmd,char *line2) {
4426 int itct,order,r,count,i,*index,er;
4427 char nm1[STRCHAR];
4428 rxnptr rxn;
4429 double rate,coeff;
4430 enum MolecState ms;
4431 listptrv vlist;
4432
4433 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
4434
4435 SCMDCHECK(line2,"missing argument");
4436 itct=strmathsscanf(line2,"%s %mlg",Varnames,Varvalues,Nvar,nm1,&coeff);
4437 SCMDCHECK(itct==2,"missing reaction name or coefficient c0");
4438 r=readrxnname(sim,nm1,&order,&rxn,&vlist,1);
4439 SCMDCHECK(r>=0,"unrecognized reaction name");
4440 line2=strnword(line2,3);
4441 rate=coeff;
4442 while(line2) {
4443 itct=strmathsscanf(line2,"%mlg %s",Varnames,Varvalues,Nvar,&coeff,nm1);
4444 SCMDCHECK(itct==2,"missing coefficient and/or species parameters");
4445 i=molstring2index1(sim,line2,&ms,&index);
4446 SCMDCHECK(i!=-1,"species is missing or cannot be read");
4447 SCMDCHECK(i!=-2,"mismatched or improper parentheses around molecule state");
4448 SCMDCHECK(i!=-3,"cannot read molecule state value");
4449 SCMDCHECK(i!=-4 || sim->ruless,"molecule name not recognized");
4450 SCMDCHECK(i!=-7,"error allocating memory");
4451 line2=strnword(line2,3);
4452 count=(i==-4)?0:molcount(sim,i,index,ms,-1);
4453 rate+=coeff*count; }
4454 if(rate<0) rate=0;
4455 if(!vlist) {
4456 er=RxnSetValue(sim,"rate",rxn,rate);
4457 SCMDCHECK(er==0,"error setting reaction rate"); }
4458 else {
4459 for(i=0;i<vlist->n;i++) {
4460 er=RxnSetValue(sim,"rate",(rxnptr) vlist->xs[i],rate);
4461 SCMDCHECK(er==0,"error setting reaction rate"); }}
4462
4463 return CMDok; }
4464
4465
4466 /* cmdexpandsystem */
cmdexpandsystem(simptr sim,cmdptr cmd,char * line2)4467 enum CMDcode cmdexpandsystem(simptr sim,cmdptr cmd,char *line2) {
4468 int itct,d,dim,s,p,c,k,i,em;
4469 double zero[DIMMAX];
4470 moleculeptr mptr;
4471 surfaceptr srf;
4472 enum PanelShape ps;
4473 enum PanelFace face;
4474 panelptr pnl;
4475 compartptr cmpt;
4476 static double expand[DIMMAX],center[DIMMAX];
4477 static int inscan=0;
4478
4479 if(inscan) goto scanportion;
4480 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
4481
4482 dim=sim->dim;
4483 SCMDCHECK(line2,"missing arguments");
4484 if(dim==1) itct=strmathsscanf(line2,"%mlg",Varnames,Varvalues,Nvar,&expand[0]);
4485 else if(dim==2) itct=strmathsscanf(line2,"%mlg %mlg",Varnames,Varvalues,Nvar,&expand[0],&expand[1]);
4486 else itct=strmathsscanf(line2,"%mlg %mlg %mlg",Varnames,Varvalues,Nvar,&expand[0],&expand[1],&expand[2]);
4487 SCMDCHECK(itct==dim,"cannot read or wrong number of expansion values");
4488 systemcenter(sim,center);
4489
4490 inscan=1;
4491 molscancmd(sim,-1,NULL,MSall,cmd,cmdexpandsystem);
4492 inscan=0;
4493
4494 if(sim->srfss) {
4495 zero[0]=zero[1]=zero[2]=0;
4496 for(s=0;s<sim->srfss->nsrf;s++) {
4497 srf=sim->srfss->srflist[s];
4498 for(ps=(enum PanelShape)0;ps<(enum PanelShape)PSMAX;ps=(enum PanelShape)(ps+1))
4499 for(p=0;p<srf->npanel[ps];p++) {
4500 pnl=srf->panels[ps][p];
4501 surftransformpanel(pnl,sim->dim,zero,center,expand); }
4502 if(srf->nemitter[PFfront] && srf->nemitter[PFback] && sim->mols)
4503 for(face=(enum PanelFace)0;face<(enum PanelFace)2;face=(enum PanelFace)(face+1))
4504 for(i=1;i<sim->mols->nspecies;i++)
4505 for(em=0;em<srf->nemitter[face][i];em++)
4506 for(d=0;d<dim;d++)
4507 srf->emitterpos[face][i][em][d]=center[d]+(srf->emitterpos[face][i][em][d]-center[d])*expand[d]; }}
4508
4509 if(sim->cmptss) {
4510 for(c=0;c<sim->cmptss->ncmpt;c++) {
4511 cmpt=sim->cmptss->cmptlist[c];
4512 for(k=0;k<cmpt->npts;k++)
4513 for(d=0;d<dim;d++)
4514 cmpt->points[k][d]=center[d]+(cmpt->points[k][d]-center[d])*expand[d]; }
4515 compartsetcondition(sim->cmptss,SCparams,0); }
4516
4517 sim->mols->touch++;
4518 return CMDok;
4519
4520 scanportion:
4521 mptr=(moleculeptr) line2;
4522 for(d=0;d<sim->dim;d++) {
4523 mptr->pos[d]=center[d]+(mptr->pos[d]-center[d])*expand[d];
4524 mptr->posx[d]=center[d]+(mptr->posx[d]-center[d])*expand[d]; }
4525 return CMDok; }
4526
4527
4528 /* cmdtranslatecmpt */
cmdtranslatecmpt(simptr sim,cmdptr cmd,char * line2)4529 enum CMDcode cmdtranslatecmpt(simptr sim,cmdptr cmd,char *line2) {
4530 int itct,dim,c,code;
4531 compartssptr cmptss;
4532 compartptr cmpt;
4533 char cname[STRCHAR];
4534 double translate[DIMMAX];
4535
4536 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
4537
4538 dim=sim->dim;
4539 cmptss=sim->cmptss;
4540 SCMDCHECK(cmptss,"no compartments defined");
4541 SCMDCHECK(line2,"first argument should be compartment name");
4542 itct=sscanf(line2,"%s",cname);
4543 SCMDCHECK(itct==1,"cannot read compartment name");
4544 c=stringfind(cmptss->cnames,cmptss->ncmpt,cname);
4545 SCMDCHECK(c>=0,"compartment name not recognized");
4546 cmpt=cmptss->cmptlist[c];
4547
4548 SCMDCHECK(line2=strnword(line2,2),"second argument should be code value");;
4549 itct=strmathsscanf(line2,"%mi",Varnames,Varvalues,Nvar,&code);
4550 SCMDCHECK(itct==1,"second argument should be code value");
4551
4552 SCMDCHECK(line2=strnword(line2,2),"missing arguments for translation amount");;
4553 if(dim==1) itct=strmathsscanf(line2,"%mlg",Varnames,Varvalues,Nvar,&translate[0]);
4554 else if(dim==2) itct=strmathsscanf(line2,"%mlg %mlg",Varnames,Varvalues,Nvar,&translate[0],&translate[1]);
4555 else itct=strmathsscanf(line2,"%mlg %mlg %mlg",Varnames,Varvalues,Nvar,&translate[0],&translate[1],&translate[2]);
4556 SCMDCHECK(itct==dim,"cannot read translation values or wrong number of them");
4557
4558 comparttranslate(sim,cmpt,code,translate);
4559
4560 return CMDok; }
4561
4562
4563 /* cmddiffusecmpt */
cmddiffusecmpt(simptr sim,cmdptr cmd,char * line2)4564 enum CMDcode cmddiffusecmpt(simptr sim,cmdptr cmd,char *line2) {
4565 int itct,dim,c,code,d,valid,pt,is,nsample,attempts;
4566 compartssptr cmptss;
4567 compartptr cmpt,cmptbound;
4568 char cname[STRCHAR];
4569 double stddev[DIMMAX],translate[DIMMAX],samplept[DIMMAX],testpt[DIMMAX],radius;
4570
4571 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
4572
4573 dim=sim->dim;
4574 cmptss=sim->cmptss;
4575 SCMDCHECK(cmptss,"no compartments defined");
4576 SCMDCHECK(line2,"first argument should be compartment name");
4577 itct=sscanf(line2,"%s",cname);
4578 SCMDCHECK(itct==1,"cannot read compartment name");
4579 c=stringfind(cmptss->cnames,cmptss->ncmpt,cname);
4580 SCMDCHECK(c>=0,"compartment name not recognized");
4581 cmpt=cmptss->cmptlist[c];
4582
4583 SCMDCHECK(line2=strnword(line2,2),"second argument should be code value");;
4584 itct=strmathsscanf(line2,"%mi",Varnames,Varvalues,Nvar,&code);
4585 SCMDCHECK(itct==1,"second argument should be code value");
4586
4587 SCMDCHECK(line2=strnword(line2,2),"missing arguments for standard deviations");;
4588 if(dim==1) itct=strmathsscanf(line2,"%mlg",Varnames,Varvalues,Nvar,&stddev[0]);
4589 else if(dim==2) itct=strmathsscanf(line2,"%mlg %mlg",Varnames,Varvalues,Nvar,&stddev[0],&stddev[1]);
4590 else itct=strmathsscanf(line2,"%mlg %mlg %mlg",Varnames,Varvalues,Nvar,&stddev[0],&stddev[1],&stddev[2]);
4591 SCMDCHECK(itct==dim,"cannot read standard deviation values or wrong number of them");
4592
4593 line2=strnword(line2,dim+1);
4594 if(line2) {
4595 itct=strmathsscanf(line2,"%s %mlg %mi",Varnames,Varvalues,Nvar,cname,&radius,&nsample);
4596 SCMDCHECK(itct==3,"cannot read bounding compartment name, radius, and/or number of samples");
4597 c=stringfind(cmptss->cnames,cmptss->ncmpt,cname);
4598 SCMDCHECK(c>=0,"bounding compartment name not recognized");
4599 cmptbound=cmptss->cmptlist[c];
4600 SCMDCHECK(cmptbound!=cmpt,"bounding compartment cannot be the same as the moving compartment");
4601 SCMDCHECK(radius>0,"bounding radius needs to be >0");
4602 SCMDCHECK(nsample>0,"number of samples needs to be >0"); }
4603 else {
4604 cmptbound=NULL;
4605 radius=0; }
4606
4607 if(!cmptbound) {
4608 for(d=0;d<dim;d++)
4609 translate[d]=stddev[d]*gaussrandD();
4610 valid=1; }
4611 else {
4612 valid=0;
4613 attempts=0;
4614 while(!valid && attempts<10) {
4615 valid=1;
4616 attempts++;
4617 for(d=0;d<dim;d++)
4618 translate[d]=stddev[d]*gaussrandD();
4619 for(pt=0;pt<cmpt->npts && valid;pt++)
4620 for(is=0;is<nsample && valid;is++) {
4621 if(dim==1) samplept[0]=(coinrandD(0.5)==0)?-radius:radius;
4622 else if(dim==2) circlerandD(samplept,radius);
4623 else sphererandCCD(samplept,radius,radius);
4624 for(d=0;d<dim;d++)
4625 testpt[d]=cmpt->points[pt][d]+samplept[d]+translate[d];
4626 if(!posincompart(sim,testpt,cmptbound,0)) {
4627 valid=0; }}}}
4628
4629 if(valid)
4630 comparttranslate(sim,cmpt,code,translate);
4631
4632 return CMDok; }
4633
4634
4635 /* cmdlongrangeforce */
cmdlongrangeforce(simptr sim,cmdptr cmd,char * line2)4636 enum CMDcode cmdlongrangeforce(simptr sim,cmdptr cmd,char *line2) {
4637 int itct,i,j,j2,ll,d,dim,wrap[DIMMAX],m,duplicate;
4638 enum MolecState ms;
4639 moleculeptr mptr,mptr2;
4640 double dt,mobility,dist,delta[DIMMAX],force;
4641 boxptr bptr;
4642
4643 static int inscan=0,i1,i2,*index1,*index2,rvar,lllo,llhi;
4644 static enum MolecState ms1,ms2,mslo,mshi;
4645 static double mobility1,mobility2,rmin,rmax,forcemag,syswidth[DIMMAX],force0[4]={0,0,0,0};
4646 static char eqstring[STRCHAR];
4647 static listptrULVD4 moleclist;
4648
4649 if(inscan) goto scanportion;
4650 if(line2 && !strcmp(line2,"cmdtype")) return CMDmanipulate;
4651
4652 i1=molstring2index1(sim,line2,&ms1,&index1);
4653 SCMDCHECK(i1!=-1,"species is missing or cannot be read");
4654 SCMDCHECK(i1!=-2,"mismatched or improper parentheses around molecule state");
4655 SCMDCHECK(i1!=-3,"cannot read molecule state value");
4656 SCMDCHECK(i1!=-4,"molecule name not recognized");
4657 SCMDCHECK(i1!=-7,"error allocating memory");
4658 line2=strnword(line2,2);
4659 i2=molstring2index1(sim,line2,&ms2,&index2);
4660 SCMDCHECK(i2!=-1,"species is missing or cannot be read");
4661 SCMDCHECK(i2!=-2,"mismatched or improper parentheses around molecule state");
4662 SCMDCHECK(i2!=-3,"cannot read molecule state value");
4663 SCMDCHECK(i2!=-4,"molecule name not recognized");
4664 SCMDCHECK(i2!=-7,"error allocating memory");
4665 line2=strnword(line2,2);
4666 SCMDCHECK(line2,"longrangeforce format: species1(state) species2(state) mobility1 mobility2 r_min r_max equation");
4667 itct=strmathsscanf(line2,"%mlg %mlg %mlg %mlg %s",Varnames,Varvalues,Nvar,&mobility1,&mobility2,&rmin,&rmax,eqstring);
4668 SCMDCHECK(itct==5,"longrangeforce format: species1(state) species2(state) mobility1 mobility2 r_min r_max equation");
4669 SCMDCHECK(rmin>0,"minimum radius needs to be >0");
4670 SCMDCHECK(rmax>=0,"maximum radius needs to be >=0");
4671
4672 dim=sim->dim;
4673 if(strhasname(eqstring,"r"))
4674 rvar=1;
4675 else {
4676 rvar=0;
4677 forcemag=strmatheval(eqstring,Varnames,Varvalues,Nvar);
4678 SCMDCHECK(forcemag==forcemag,"cannot compute equation value"); }
4679
4680 lllo=llhi=-1;
4681 if(ms2<MSMAX) {
4682 mslo=ms2;
4683 mshi=(enum MolecState)(ms2+1); }
4684 else {
4685 mslo=(enum MolecState) 0;
4686 mshi=(enum MolecState) MSMAX; }
4687 for(i=0;i<index2[PDnresults];i++)
4688 for(ms=mslo;ms<mshi;ms=(enum MolecState)(ms+1)) {
4689 ll=sim->mols->listlookup[index2[PDMAX+i]][ms];
4690 if(ll<lllo || lllo==-1) lllo=ll;
4691 if(ll>=llhi || llhi==-1) llhi=ll+1; }
4692
4693 if(!cmd->v1) {
4694 cmd->v1=(void*) ListAllocULVD4(256); // start with 256 molecules, and expand as needed
4695 SCMDCHECK(cmd->v1,"memory allocation error");
4696 cmd->freefn=&cmdListULVD4free; }
4697
4698 moleclist=(listptrULVD4) cmd->v1;
4699 for(j=0;j<moleclist->n;j++) {
4700 moleclist->datad4[j][3]=0;
4701 for(d=0;d<dim;d++)
4702 moleclist->datad4[j][d]=0; }
4703
4704 for(d=0;d<dim;d++)
4705 syswidth[d]=sim->wlist[2*d+1]->pos-sim->wlist[2*d]->pos;
4706
4707 inscan=1;
4708 molscancmd(sim,i1,index1,ms1,cmd,cmdlongrangeforce);
4709 inscan=0;
4710
4711 dt=sim->dt;
4712 for(i=0;i<moleclist->n;i++) {
4713 if(moleclist->datad4[i][3]==0)
4714 moleclist->datav[i]=NULL;
4715 else {
4716 mptr=(moleculeptr) moleclist->datav[i];
4717 mobility=molismatch(mptr,i1,index1,ms1)?mobility1:mobility2;
4718 for(d=0;d<dim;d++)
4719 mptr->pos[d]+=moleclist->datad4[i][d]*mobility*dt; }}
4720 List_CleanULVD4(moleclist);
4721
4722 return CMDok;
4723
4724 scanportion:
4725 dim=sim->dim;
4726 mptr=(moleculeptr) line2;
4727 j=ListInsertItemULVD4(moleclist,mptr->serno,(void*)mptr,force0,1);
4728 SCMDCHECK(j>=0,"failed to allocate memory");
4729 moleclist->datad4[j][3]=1;
4730 bptr=boxscansphere(sim,mptr->pos,rmax,NULL,wrap);
4731 while(bptr) {
4732 for(ll=lllo;ll<llhi;ll++)
4733 for(m=0;m<bptr->nmol[ll];m++) {
4734 mptr2=bptr->mol[ll][m];
4735 if(molismatch(mptr2,i2,index2,ms2) && mptr2!=mptr) {
4736 dist=0;
4737 for(d=0;d<dim;d++) {
4738 delta[d]=mptr2->pos[d]+wrap[d]*syswidth[d]-mptr->pos[d];
4739 dist+=delta[d]*delta[d]; }
4740 dist=sqrt(dist);
4741 if(dist>=rmin && dist<=rmax) {
4742 if(rvar) {
4743 simsetvariable(sim,"r",dist);
4744 forcemag=strmatheval(eqstring,Varnames,Varvalues,Nvar); }
4745 duplicate=molismatch(mptr2,i1,index1,ms1);
4746 if(!duplicate) {
4747 j2=ListInsertItemULVD4(moleclist,mptr2->serno,(void*)mptr2,force0,1);
4748 moleclist->datad4[j2][3]=1;
4749 j=ListInsertItemULVD4(moleclist,mptr->serno,NULL,NULL,0); } // in case adding j2 moved the position for j
4750 for(d=0;d<dim;d++) {
4751 force=delta[d]/dist*forcemag;
4752 moleclist->datad4[j][d]-=force;
4753 if(!duplicate)
4754 moleclist->datad4[j2][d]+=force; }}}}
4755 bptr=boxscansphere(sim,mptr->pos,rmax,bptr,wrap); }
4756 return CMDok; }
4757
4758
4759
4760
4761
4762 /**********************************************************/
4763 /******************* internal routines ********************/
4764 /**********************************************************/
4765
4766
4767 /* cmdv1free */
cmdv1free(cmdptr cmd)4768 void cmdv1free(cmdptr cmd) {
4769 free(cmd->v1);
4770 cmd->v1=NULL;
4771 return; }
4772
4773
4774 /* cmdv1v2free */
cmdv1v2free(cmdptr cmd)4775 void cmdv1v2free(cmdptr cmd) {
4776 free(cmd->v1);
4777 free(cmd->v2);
4778 cmd->v1=NULL;
4779 cmd->v2=NULL;
4780 return; }
4781
4782
4783 /* cmdListULVD3free */
cmdListULVD4free(cmdptr cmd)4784 void cmdListULVD4free(cmdptr cmd) {
4785 ListFreeULVD4((listptrULVD4) cmd->v1);
4786 cmd->v1=NULL;
4787 return; }
4788
4789
4790 /* conditionalcmdtype */
conditionalcmdtype(simptr sim,cmdptr cmd,int nparam)4791 enum CMDcode conditionalcmdtype(simptr sim,cmdptr cmd,int nparam) {
4792 char string[STRCHAR],*strptr;
4793 enum CMDcode ans;
4794
4795 if(!cmd->str) return CMDnone;
4796 strptr=strnword(cmd->str,nparam+2);
4797 if(!strptr) return CMDnone;
4798 strcpy(string,strptr);
4799 strptr=cmd->str;
4800 cmd->str=string;
4801 ans=scmdcmdtype(sim->cmds,cmd);
4802 cmd->str=strptr;
4803 return ans; }
4804
4805
4806 /* insideecoli */
insideecoli(double * pos,double * ofst,double rad,double length)4807 int insideecoli(double *pos,double *ofst,double rad,double length) {
4808 double dist;
4809
4810 dist=(pos[1]-ofst[1])*(pos[1]-ofst[1])+(pos[2]-ofst[2])*(pos[2]-ofst[2]);
4811 if(pos[0]-ofst[0]<rad) dist+=(pos[0]-ofst[0]-rad)*(pos[0]-ofst[0]-rad);
4812 else if(pos[0]-ofst[0]>length-rad) dist+=(pos[0]-ofst[0]-length+rad)*(pos[0]-ofst[0]-length+rad);
4813 return dist<rad*rad; }
4814
4815
4816 /* putinecoli */
putinecoli(double * pos,double * ofst,double rad,double length)4817 void putinecoli(double *pos,double *ofst,double rad,double length) {
4818 double dist;
4819
4820 dist=(pos[1]-ofst[1])*(pos[1]-ofst[1])+(pos[2]-ofst[2])*(pos[2]-ofst[2]);
4821 if(pos[0]-ofst[0]<rad) {
4822 dist+=(pos[0]-ofst[0]-rad)*(pos[0]-ofst[0]-rad);
4823 dist=sqrt(rad*rad/dist);
4824 pos[0]=ofst[0]+rad+dist*(pos[0]-ofst[0]-rad); }
4825 else if(pos[0]-ofst[0]>length-rad) {
4826 dist+=(pos[0]-ofst[0]-length+rad)*(pos[0]-ofst[0]-length+rad);
4827 dist=sqrt(rad*rad/dist);
4828 pos[0]=ofst[0]+length-rad+dist*(pos[0]-ofst[0]-length+rad); }
4829 else
4830 dist=sqrt(rad*rad/dist);
4831 pos[1]=ofst[1]+dist*(pos[1]-ofst[1]);
4832 pos[2]=ofst[2]+dist*(pos[2]-ofst[2]);
4833 return; }
4834
4835
4836 /* molinpanels */
molinpanels(simptr sim,moleculeptr mptr,int s,enum PanelShape ps)4837 int molinpanels(simptr sim,moleculeptr mptr,int s,enum PanelShape ps) {
4838 int p,dim,npnl;
4839 double *pos;
4840 panelptr *pnls,pnl;
4841
4842 if(ps!=PSsph) return 0;
4843
4844 if(s<0) {
4845 for(s=0;s<sim->srfss->nsrf;s++)
4846 if(molinpanels(sim,mptr,s,ps)) return 1;
4847 return 0; }
4848
4849 dim=sim->dim;
4850 pnls=sim->srfss->srflist[s]->panels[ps];
4851 npnl=sim->srfss->srflist[s]->npanel[ps];
4852 pos=mptr->pos;
4853 for(p=0;p<npnl;p++) {
4854 pnl=pnls[p];
4855 if(Geo_PtInSphere(pos,pnl->point[0],pnl->point[1][0],dim)) return 1; }
4856 return 0; }
4857
4858