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 <float.h>
9 #include <stdio.h>
10 #include <string.h>
11 #include "Geometry.h"
12 #include "math2.h"
13 #include "random2.h"
14 #include "smoldyn.h"
15 #include "smoldynfuncs.h"
16 #include "Zn.h"
17
18 /******************************************************************************/
19 /******************************** Virtual boxes *******************************/
20 /******************************************************************************/
21
22
23 /******************************************************************************/
24 /****************************** Local declarations ****************************/
25 /******************************************************************************/
26
27 // low level utilities
28 int panelinbox(simptr sim,panelptr pnl,boxptr bptr);
29
30 // memory management
31 boxptr boxalloc(int dim,int nlist);
32 int expandbox(boxptr bptr,int n,int ll);
33 int expandboxpanels(boxptr bptr,int n);
34 void boxfree(boxptr bptr,int nlist);
35 boxptr *boxesalloc(int nbox,int dim,int nlist);
36 void boxesfree(boxptr *blist,int nbox,int nlist);
37 boxssptr boxssalloc(int dim);
38
39 // data structure output
40
41 // structure set up
42 int boxesupdateparams(simptr sim);
43 int boxesupdatelists(simptr sim);
44
45 // core simulation functions
46
47
48 /******************************************************************************/
49 /***************************** low level utilities ****************************/
50 /******************************************************************************/
51
52
53 /* box2pos */
box2pos(simptr sim,boxptr bptr,double * poslo,double * poshi)54 void box2pos(simptr sim,boxptr bptr,double *poslo,double *poshi) {
55 int d,dim;
56 double *size,*min;
57
58 dim=sim->dim;
59 size=sim->boxs->size;
60 min=sim->boxs->min;
61 if(poslo) for(d=0;d<dim;d++) poslo[d]=min[d]+bptr->indx[d]*size[d];
62 if(poshi) for(d=0;d<dim;d++) poshi[d]=min[d]+(bptr->indx[d]+1)*size[d];
63 return; }
64
65
66 /* pos2box */
pos2box(simptr sim,const double * pos)67 boxptr pos2box(simptr sim,const double *pos) {
68 int b,d,indx,dim;
69 boxssptr boxs;
70
71 dim=sim->dim;
72 boxs=sim->boxs;
73 b=0;
74 for(d=0;d<dim;d++) {
75 indx=(int)((pos[d]-boxs->min[d])/boxs->size[d]);
76 if(indx<0) indx=0;
77 else if(indx>=boxs->side[d]) indx=boxs->side[d]-1;
78 b=boxs->side[d]*b+indx; }
79 return boxs->blist[b]; }
80
81
82 /* boxrandpos */
boxrandpos(simptr sim,double * pos,boxptr bptr)83 void boxrandpos(simptr sim,double *pos,boxptr bptr) {
84 int d;
85 double *min,*size;
86
87 min=sim->boxs->min;
88 size=sim->boxs->size;
89 for(d=0;d<sim->dim;d++)
90 pos[d]=unirandCCD(min[d]+size[d]*bptr->indx[d],min[d]+size[d]*(bptr->indx[d]+1));
91 return; }
92
93
94 /* panelinbox */
panelinbox(simptr sim,panelptr pnl,boxptr bptr)95 int panelinbox(simptr sim,panelptr pnl,boxptr bptr) {
96 int dim,d,cross;
97 double v1[DIMMAX],v2[DIMMAX],v3[DIMMAX],v4[DIMMAX],**point,*front;
98
99 dim=sim->dim;
100 box2pos(sim,bptr,v1,v2); // v1 and v2 are set to corners of box
101 for(d=0;d<dim;d++) {
102 if(bptr->indx[d]==0) v1[d]=-VERYLARGE;
103 if(bptr->indx[d]==sim->boxs->side[d]-1) v2[d]=VERYLARGE; }
104 point=pnl->point;
105 front=pnl->front;
106
107 if(pnl->ps==PSrect) {
108 if(dim==1) cross=Geo_PtInSlab(v1,v2,point[0],dim);
109 else if(dim==2) {
110 v3[0]=front[1]==0?1:0;
111 v3[1]=front[1]==0?0:1;
112 cross=Geo_LineXaabb2(point[0],point[1],v3,v1,v2); }
113 else cross=Geo_RectXaabb3(point[0],point[1],point[3],point[0],v1,v2); }
114 else if(pnl->ps==PStri) {
115 if(dim==1) cross=Geo_PtInSlab(v1,v2,point[0],dim);
116 else if(dim==2) cross=Geo_LineXaabb2(point[0],point[1],front,v1,v2);
117 else cross=Geo_TriXaabb3(point[0],point[1],point[2],front,v1,v2); }
118 else if(pnl->ps==PSsph) {
119 if(dim==1) {
120 if((point[0][0]-point[1][0]<v1[0] || point[0][0]-point[1][0]>=v2[0]) && (point[0][0]+point[1][0]<v1[0] || point[0][0]+point[1][0]>=v2[0])) cross=0;
121 else cross=1; }
122 else if(dim==2) cross=Geo_CircleXaabb2(point[0],point[1][0],v1,v2);
123 else cross=Geo_SphsXaabb3(point[0],point[1][0],v1,v2); }
124 else if(pnl->ps==PScyl) {
125 if(dim==2) {
126 v3[0]=point[0][0]+point[2][0]*front[0];
127 v3[1]=point[0][1]+point[2][0]*front[1];
128 v4[0]=point[1][0]+point[2][0]*front[0];
129 v4[1]=point[1][1]+point[2][0]*front[1];
130 cross=Geo_LineXaabb2(v3,v4,front,v1,v2);
131 if(!cross) {
132 v3[0]=point[0][0]-point[2][0]*front[0];
133 v3[1]=point[0][1]-point[2][0]*front[1];
134 v4[0]=point[1][0]-point[2][0]*front[0];
135 v4[1]=point[1][1]-point[2][0]*front[1];
136 cross=Geo_LineXaabb2(v3,v4,front,v1,v2); }}
137 else {
138 cross=Geo_CylsXaabb3(point[0],point[1],point[2][0],v1,v2); }}
139 else if(pnl->ps==PShemi) {
140 if(dim==2) cross=Geo_SemicXaabb2(point[0],point[1][0],point[2],v1,v2);
141 else cross=Geo_HemisXaabb3(point[0],point[1][0],point[2],v1,v2); }
142 else if(pnl->ps==PSdisk) {
143 if(dim==2) {
144 v3[0]=point[0][0]+point[1][0]*front[1];
145 v3[1]=point[0][1]-point[1][0]*front[0];
146 v4[0]=point[0][0]-point[1][0]*front[1];
147 v4[1]=point[0][1]+point[1][0]*front[0];
148 cross=Geo_LineXaabb2(v3,v4,front,v1,v2); }
149 else cross=Geo_DiskXaabb3(point[0],point[1][0],front,v1,v2); }
150 else cross=0;
151
152 return cross; }
153
154
155 /* boxaddmol */
boxaddmol(moleculeptr mptr,int ll)156 int boxaddmol(moleculeptr mptr,int ll) {
157 boxptr bptr;
158
159 bptr=mptr->box;
160 if(bptr->nmol[ll]==bptr->maxmol[ll])
161 if(expandbox(bptr,bptr->maxmol[ll]+1,ll)) return 1;
162 bptr->mol[ll][bptr->nmol[ll]++]=mptr;
163 return 0; }
164
165
166 /* boxremovemol */
boxremovemol(moleculeptr mptr,int ll)167 void boxremovemol(moleculeptr mptr,int ll) {
168 int m;
169 boxptr bptr;
170
171 bptr=mptr->box;
172 for(m=bptr->nmol[ll]-1;m>=0 && bptr->mol[ll][m]!=mptr;m--);
173 if(m>=0) bptr->mol[ll][m]=bptr->mol[ll][--bptr->nmol[ll]];
174 mptr->box=NULL;
175 return; }
176
177
178 /* boxscansphere */
boxscansphere(simptr sim,const double * pos,double radius,boxptr bptr,int * wrap)179 boxptr boxscansphere(simptr sim,const double *pos,double radius,boxptr bptr,int *wrap) {
180 boxssptr boxs;
181 double boxmin[DIMMAX],boxmax[DIMMAX],v1[DIMMAX],dist;
182 int dim,d,b,index[DIMMAX],done,diam,keepgoing;
183 static int boxdiameter,startindex[DIMMAX],deltaindex[DIMMAX];
184
185 boxs=sim->boxs;
186 dim=sim->dim;
187
188 if(bptr==NULL) { // first call
189 boxdiameter=0;
190 for(d=0;d<dim;d++) {
191 startindex[d]=(int) floor(((pos[d]-boxs->min[d])-radius)/boxs->size[d]); // first box to consider
192 diam=(int) ceil(((pos[d]-boxs->min[d])+radius)/boxs->size[d])-startindex[d];
193 if(diam>boxdiameter) boxdiameter=diam; // largest diameter in boxes, giving count range
194 deltaindex[d]=0; }} // working index relative to startindex
195 else { // subsequent calls
196 done=Zn_incrementcounter(deltaindex,dim,boxdiameter);
197 if(done) return NULL; }
198
199 keepgoing=1;
200 while(keepgoing) {
201 keepgoing=0; // expect a good box
202 for(d=0;d<dim;d++) {
203 index[d]=startindex[d]+deltaindex[d]; // index of box, ignoring edges and wrapping
204 boxmin[d]=boxs->min[d]+boxs->size[d]*index[d]; // min corner of box, ignoring edges and wrapping
205 boxmax[d]=boxs->min[d]+boxs->size[d]*(index[d]+1); // max corner of box, ignoring edges and wrapping
206 wrap[d]=0;
207 if(index[d]<0) {
208 if(sim->wlist[2*d]->type=='p')
209 while(index[d]<0) {
210 index[d]+=boxs->side[d];
211 wrap[d]--; }
212 else
213 keepgoing=1; }
214 else if(index[d]>=boxs->side[d]) {
215 if(sim->wlist[2*d+1]->type=='p')
216 while(index[d]>=boxs->side[d]) {
217 index[d]-=boxs->side[d];
218 wrap[d]++; }
219 else
220 keepgoing=1; }}
221 dist=Geo_NearestAabbPt(boxmin,boxmax,dim,pos,v1);
222 if(dist>radius)
223 keepgoing=1;
224 if(keepgoing) {
225 done=Zn_incrementcounter(deltaindex,dim,boxdiameter);
226 if(done) return NULL; }}
227
228 b=0; // compute the box number with this index
229 for(d=0;d<dim;d++)
230 b=boxs->side[d]*b+index[d];
231 bptr=boxs->blist[b];
232
233 return bptr; }
234
235
236 /* boxdebug */
boxdebug(simptr sim)237 int boxdebug(simptr sim) {
238 int er,ll,m,mb,b;
239 molssptr mols;
240 boxssptr boxs;
241 moleculeptr mptr;
242 boxptr bptr;
243 char string[STRCHAR];
244
245 er=0;
246 mols=sim->mols;
247 boxs=sim->boxs;
248
249 for(b=0;b<boxs->nbox;b++) {
250 bptr=boxs->blist[b];
251 for(ll=0;ll<mols->nlist;ll++) {
252 printf("Box %p list %i:",bptr,ll);
253 for(mb=0;mb<bptr->nmol[ll];mb++)
254 printf(" %s",molserno2string(bptr->mol[ll][mb]->serno,string));
255 printf("\n"); }}
256
257
258 for(ll=0;ll<mols->nlist;ll++) // check that molecules are in the boxes they think they are in
259 for(m=0;m<mols->nl[ll];m++) {
260 mptr=mols->live[ll][m];
261 bptr=mptr->box;
262 if(!bptr) {er++;printf("BUG: molecule %s has box value set to NULL\n",molserno2string(mptr->serno,string)); }
263 else {
264 for(mb=0;mb<bptr->nmol[ll] && bptr->mol[ll][mb]!=mptr;mb++);
265 if(mb==bptr->nmol[ll]) {er++;printf("BUG: molecule %s thinks it's in box %p but isn't\n",molserno2string(mptr->serno,string),bptr); } }}
266
267 for(b=0;b<boxs->nbox;b++) {
268 bptr=boxs->blist[b];
269 for(ll=0;ll<mols->nlist;ll++)
270 for(mb=0;mb<bptr->nmol[ll];mb++) {
271 mptr=bptr->mol[ll][mb];
272 if(mptr->box!=bptr) {er++;printf("BUG: molecule %s thinks it's in box %p but is really in box %p\n",molserno2string(mptr->serno,string),mptr->box,bptr);} }}
273
274 return er; }
275
276
277 /******************************************************************************/
278 /****************************** memory management *****************************/
279 /******************************************************************************/
280
281 /* boxalloc */
boxalloc(int dim,int nlist)282 boxptr boxalloc(int dim,int nlist) {
283 boxptr bptr;
284 int d,ll;
285
286 bptr=NULL;
287 CHECKMEM(bptr=(boxptr) malloc(sizeof(struct boxstruct)));
288 bptr->indx=NULL;
289 bptr->nneigh=0;
290 bptr->midneigh=0;
291 bptr->neigh=NULL;
292 bptr->wpneigh=NULL;
293 bptr->nwall=0;
294 bptr->wlist=NULL;
295 bptr->maxpanel=0;
296 bptr->npanel=0;
297 bptr->panel=NULL;
298 bptr->maxmol=NULL;
299 bptr->nmol=NULL;
300 bptr->mol=NULL;
301
302 CHECKMEM(bptr->indx=(int*) calloc(dim,sizeof(int)));
303 for(d=0;d<dim;d++) bptr->indx[d]=0;
304 if(nlist) {
305 CHECKMEM(bptr->maxmol=(int*) calloc(nlist,sizeof(int)));
306 for(ll=0;ll<nlist;ll++) bptr->maxmol[ll]=0;
307 CHECKMEM(bptr->nmol=(int*) calloc(nlist,sizeof(int)));
308 for(ll=0;ll<nlist;ll++) bptr->nmol[ll]=0;
309 CHECKMEM(bptr->mol=(moleculeptr**) calloc(nlist,sizeof(moleculeptr*)));
310 for(ll=0;ll<nlist;ll++) bptr->mol[ll]=NULL; }
311
312 return bptr;
313
314 failure:
315 boxfree(bptr,nlist);
316 simLog(NULL,10,"Failed to allocate memory in boxalloc");
317 return NULL; }
318
319
320 /* expandbox */
expandbox(boxptr bptr,int n,int ll)321 int expandbox(boxptr bptr,int n,int ll) {
322 moleculeptr *mlist;
323 int m,maxmol,mn;
324
325 maxmol=bptr->maxmol[ll]+n;
326 if(maxmol>0) {
327 mlist=(moleculeptr*) calloc(maxmol,sizeof(moleculeptr));
328 if(!mlist) return 1;
329 mn=(n>0)?bptr->maxmol[ll]:maxmol;
330 for(m=0;m<mn;m++) mlist[m]=bptr->mol[ll][m]; }
331 else {
332 maxmol=0;
333 mlist=NULL; }
334 free(bptr->mol[ll]);
335 bptr->mol[ll]=mlist;
336 bptr->maxmol[ll]=maxmol;
337 if(bptr->nmol[ll]>maxmol) bptr->nmol[ll]=maxmol;
338 return 0;}
339
340
341 /* expandboxpanels */
expandboxpanels(boxptr bptr,int n)342 int expandboxpanels(boxptr bptr,int n) {
343 int maxpanel,p;
344 panelptr *panel;
345
346 if(n<=0) return 0;
347 maxpanel=n+bptr->maxpanel;
348 panel=(panelptr*) calloc(maxpanel,sizeof(panelptr));
349 if(!panel) return 1;
350 for(p=0;p<bptr->npanel;p++)
351 panel[p]=bptr->panel[p];
352 for(;p<maxpanel;p++)
353 panel[p]=NULL;
354 free(bptr->panel);
355 bptr->panel=panel;
356 bptr->maxpanel=maxpanel;
357 return 0; }
358
359
360 /* boxfree */
boxfree(boxptr bptr,int nlist)361 void boxfree(boxptr bptr,int nlist) {
362 int ll;
363
364 if(!bptr) return;
365 if(bptr->mol) {
366 for(ll=0;ll<nlist;ll++)
367 free(bptr->mol[ll]); }
368 free(bptr->mol);
369 free(bptr->nmol);
370 free(bptr->maxmol);
371 free(bptr->panel);
372 free(bptr->wlist);
373 free(bptr->wpneigh);
374 free(bptr->neigh);
375 free(bptr->indx);
376 free(bptr);
377 return; }
378
379
380 /* boxesalloc */
boxesalloc(int nbox,int dim,int nlist)381 boxptr *boxesalloc(int nbox,int dim,int nlist) {
382 int b;
383 boxptr *blist;
384
385 blist=NULL;
386 CHECKMEM(blist=(boxptr*) calloc(nbox,sizeof(boxptr)));
387 for(b=0;b<nbox;b++) blist[b]=NULL;
388 for(b=0;b<nbox;b++)
389 CHECKMEM(blist[b]=boxalloc(dim,nlist));
390 return blist;
391
392 failure:
393 boxesfree(blist,nbox,nlist);
394 simLog(NULL,10,"Failed to allocate memory in boxesalloc");
395 return NULL; }
396
397
398 /* boxesfree */
boxesfree(boxptr * blist,int nbox,int nlist)399 void boxesfree(boxptr *blist,int nbox,int nlist) {
400 int b;
401
402 if(!blist) return;
403 for(b=0;b<nbox;b++) boxfree(blist[b],nlist);
404 free(blist);
405 return; }
406
407
408 /* boxssalloc */
boxssalloc(int dim)409 boxssptr boxssalloc(int dim) {
410 boxssptr boxs;
411 int d;
412
413 boxs=NULL;
414 CHECKMEM(boxs=(boxssptr) malloc(sizeof(struct boxsuperstruct)));
415 boxs->condition=SCinit;
416 boxs->sim=NULL;
417 boxs->nlist=0;
418 boxs->mpbox=0;
419 boxs->boxsize=0;
420 boxs->boxvol=0;
421 boxs->nbox=0;
422 boxs->side=NULL;
423 boxs->min=NULL;
424 boxs->size=NULL;
425 boxs->blist=NULL;
426
427 CHECKMEM(boxs->side=(int*) calloc(dim,sizeof(int)));
428 for(d=0;d<dim;d++) boxs->side[d]=0;
429 CHECKMEM(boxs->min=(double*) calloc(dim,sizeof(double)));
430 for(d=0;d<dim;d++) boxs->min[d]=0;
431 CHECKMEM(boxs->size=(double*) calloc(dim,sizeof(double)));
432 for(d=0;d<dim;d++) boxs->size[d]=0;
433 return boxs;
434
435 failure:
436 boxssfree(boxs);
437 simLog(NULL,10,"Failed to allocate memory in boxssalloc");
438 return NULL; }
439
440
441 /* boxssfree */
boxssfree(boxssptr boxs)442 void boxssfree(boxssptr boxs) {
443 if(!boxs) return;
444 boxesfree(boxs->blist,boxs->nbox,boxs->nlist);
445 free(boxs->size);
446 free(boxs->min);
447 free(boxs->side);
448 free(boxs);
449 return; }
450
451
452 /******************************************************************************/
453 /*************************** data structure output ****************************/
454 /******************************************************************************/
455
456
457 /* boxoutput */
boxoutput(boxssptr boxs,int blo,int bhi,int dim)458 void boxoutput(boxssptr boxs,int blo,int bhi,int dim) {
459 int b,b2,b3,p,d,ll;
460 boxptr bptr;
461 simptr sim;
462
463 sim=boxs->sim;
464 simLog(sim,2,"INDIVIDUAL BOX PARAMETERS\n");
465 if(!boxs) {
466 simLog(sim,2," No box superstructure defined\n\n");
467 return; }
468 if(bhi<0 || bhi>boxs->nbox) bhi=boxs->nbox;
469 for(b=blo;b<bhi;b++) {
470 bptr=boxs->blist[b];
471 simLog(sim,2," Box %i: indx=(",b);
472 for(d=0;d<dim-1;d++) simLog(sim,2,"%i,",bptr->indx[d]);
473 simLog(sim,2,"%i), nwall=%i\n",bptr->indx[d],bptr->nwall);
474
475 simLog(sim,2," nneigh=%i midneigh=%i\n",bptr->nneigh,bptr->midneigh);
476 if(bptr->neigh) {
477 simLog(sim,2," neighbors:");
478 for(b2=0;b2<bptr->nneigh;b2++) {
479 b3=indx2addZV(bptr->neigh[b2]->indx,boxs->side,dim);
480 simLog(sim,2," %i",b3); }
481 simLog(sim,2,"\n"); }
482 if(bptr->wpneigh) {
483 simLog(sim,2," wrap code:");
484 for(b2=0;b2<bptr->nneigh;b2++) simLog(sim,2," %i",bptr->wpneigh[b2]);
485 simLog(sim,2,"\n"); }
486
487 simLog(sim,2," %i panels",bptr->npanel);
488 if(bptr->npanel) {
489 simLog(sim,2,": ");
490 for(p=0;p<bptr->npanel;p++) {
491 simLog(sim,2," %s",bptr->panel[p]->pname); }}
492 simLog(sim,2,"\n");
493
494 simLog(sim,2," %i live lists:\n",boxs->nlist);
495 simLog(sim,2," max:");
496 for(ll=0;ll<boxs->nlist;ll++) simLog(sim,2," %i",bptr->maxmol[ll]);
497 simLog(sim,2,"\n size:");
498 for(ll=0;ll<boxs->nlist;ll++) simLog(sim,2," %i",bptr->nmol[ll]);
499 simLog(sim,2,"\n"); }
500
501 if(b<boxs->nbox) simLog(sim,2," ...\n");
502 simLog(sim,2,"\n");
503 return; }
504
505
506 /* boxssoutput */
boxssoutput(simptr sim)507 void boxssoutput(simptr sim) {
508 int dim,d,ll;
509 boxssptr boxs;
510 double flt1;
511
512 simLog(sim,2,"VIRTUAL BOX PARAMETERS\n");
513 if(!sim || !sim->boxs) {
514 simLog(sim,2," No box superstructure defined\n\n");
515 return; }
516 dim=sim->dim;
517 boxs=sim->boxs;
518 simLog(sim,2," %i boxes\n",boxs->nbox);
519 simLog(sim,2," Number of boxes on each side:");
520 for(d=0;d<dim;d++) simLog(sim,2," %i",boxs->side[d]);
521 simLog(sim,2,"\n");
522 simLog(sim,1," Minimum box position: ");
523 for(d=0;d<dim;d++) simLog(sim,1," %g",boxs->min[d]);
524 simLog(sim,1,"\n");
525 if(boxs->boxsize) simLog(sim,2," Requested box width: %g\n",boxs->boxsize);
526 if(boxs->mpbox) simLog(sim,2," Requested molecules per box: %g\n",boxs->mpbox);
527 simLog(sim,2," Box dimensions: ");
528 for(d=0;d<dim;d++) simLog(sim,2," %g",boxs->size[d]);
529 simLog(sim,2,"\n");
530 if(boxs->boxvol>0)
531 simLog(sim,2," Box volumes: %g\n",boxs->boxvol);
532 else
533 simLog(sim,2," Box volumes not computed\n");
534
535 if(sim->mols) {
536 flt1=0;
537 for(ll=0;ll<sim->mols->nlist;ll++)
538 if(sim->mols->listtype[ll]==MLTsystem) flt1+=sim->mols->nl[ll];
539 flt1/=boxs->nbox;
540 simLog(sim,2," Molecules per box= %g\n",flt1);
541 simLog(sim,2,"\n"); }
542
543 return; }
544
545
546 /* checkboxparams */
checkboxparams(simptr sim,int * warnptr)547 int checkboxparams(simptr sim,int *warnptr) {
548 int error,warn,b,dim,nmolec,ll;
549 boxssptr boxs;
550 double mpbox;
551 char string[STRCHAR];
552 boxptr bptr;
553
554 error=warn=0;
555 boxs=sim->boxs;
556 dim=sim->dim;
557 if(!boxs) {
558 warn++;
559 simLog(sim,9," WARNING: no box structure defined\n");
560 if(warnptr) *warnptr=warn;
561 return 0; }
562
563 if(boxs->condition!=SCok) {
564 warn++;
565 simLog(sim,7," WARNING: box structure %s\n",simsc2string(boxs->condition,string)); }
566
567 if(boxs->mpbox>100) { // mpbox value
568 warn++;
569 simLog(sim,5," WARNING: requested molecules per box, %g, is very high\n",boxs->mpbox); }
570 else if(boxs->mpbox>0 && boxs->mpbox<1) {
571 warn++;
572 simLog(sim,5," WARNING: requested molecules per box, %g, is very low\n",boxs->mpbox); }
573 mpbox=boxs->mpbox>0?boxs->mpbox:10;
574
575 for(b=0;b<boxs->nbox;b++) {
576 bptr=boxs->blist[b];
577 if(sim->mols) {
578 nmolec=0; // actual molecs per box
579 for(ll=0;ll<sim->mols->nlist;ll++) nmolec+=bptr->nmol[ll];
580 if(nmolec>10*mpbox) {
581 warn++;
582 simLog(sim,5," WARNING: box (%s) has %i molecules in it, which is very high\n",Zn_vect2csvstring(bptr->indx,dim,string),nmolec); }}
583
584 if(bptr->npanel>20) {
585 warn++;
586 simLog(sim,5," WARNING: box (%s) has %i panels in it, which is very high\n",Zn_vect2csvstring(bptr->indx,dim,string),bptr->npanel); }}
587
588 if(warnptr) *warnptr=warn;
589 return error; }
590
591
592
593 /******************************************************************************/
594 /********************************* structure set up ***************************/
595 /******************************************************************************/
596
597
598 /* boxsetcondition */
boxsetcondition(boxssptr boxs,enum StructCond cond,int upgrade)599 void boxsetcondition(boxssptr boxs,enum StructCond cond,int upgrade) {
600 if(!boxs) return;
601 if(upgrade==0 && boxs->condition>cond) boxs->condition=cond;
602 else if(upgrade==1 && boxs->condition<cond) boxs->condition=cond;
603 else if(upgrade==2) boxs->condition=cond;
604 if(boxs->sim && boxs->condition<boxs->sim->condition) {
605 cond=boxs->condition;
606 simsetcondition(boxs->sim,cond==SCinit?SClists:cond,0); }
607 return; }
608
609
610 /* boxsetsize */
boxsetsize(simptr sim,const char * info,double val)611 int boxsetsize(simptr sim,const char *info,double val) {
612 boxssptr boxs;
613
614 if(val<=0) return 2;
615 if(!sim->boxs) {
616 if(!sim->dim) return 3;
617 boxs=boxssalloc(sim->dim);
618 if(!boxs) return 1;
619 boxs->sim=sim;
620 sim->boxs=boxs;
621 boxsetcondition(boxs,SCinit,0); }
622 else
623 boxs=sim->boxs;
624 if(!strcmp(info,"molperbox")) boxs->mpbox=val;
625 else if(!strcmp(info,"boxsize")) boxs->boxsize=val;
626 else return 2;
627 boxsetcondition(boxs,SClists,0);
628 return 0; }
629
630
631 /* boxesupdateparams */
boxesupdateparams(simptr sim)632 int boxesupdateparams(simptr sim) {
633 int m,mlo,mhi,nbox,b,ll,ll1,mxml,er,npanel;
634 boxssptr boxs;
635 boxptr *blist,bptr;
636 int nsrf,s,p;
637 surfaceptr srf;
638 moleculeptr mptr,*mlist;
639 enum PanelShape ps;
640 panelptr pnl;
641
642 boxs=sim->boxs;
643 nbox=boxs->nbox;
644 blist=boxs->blist;
645
646 if(sim->srfss) {
647 for(b=0;b<nbox;b++)
648 blist[b]->npanel=0;
649 nsrf=sim->srfss->nsrf;
650 for(b=0;b<nbox;b++) { // box->npanel, panel
651 bptr=blist[b];
652 npanel=0;
653 for(s=0;s<nsrf;s++) {
654 srf=sim->srfss->srflist[s];
655 for(ps=(enum PanelShape)0;ps<PSMAX;ps=(enum PanelShape)(ps+1))
656 for(p=0;p<srf->npanel[ps];p++) {
657 pnl=srf->panels[ps][p];
658 if(panelinbox(sim,pnl,bptr)) npanel++; }}
659 if(npanel && npanel>bptr->maxpanel) {
660 er=expandboxpanels(bptr,npanel-bptr->maxpanel);
661 if(er) return 1; }
662 if(npanel) {
663 for(s=0;s<nsrf;s++) {
664 srf=sim->srfss->srflist[s];
665 for(ps=(enum PanelShape)0;ps<PSMAX;ps=(enum PanelShape)(ps+1))
666 for(p=0;p<srf->npanel[ps];p++)
667 if(panelinbox(sim,srf->panels[ps][p],bptr))
668 bptr->panel[bptr->npanel++]=srf->panels[ps][p]; }}}}
669
670 if(sim->mols) { // mptr->box, box->maxmol, nmol, mol
671 if(sim->mols->condition<SCparams) return 2;
672 for(b=0;b<nbox;b++) // clear out molecule lists in boxes
673 for(ll=0;ll<boxs->nlist;ll++)
674 blist[b]->nmol[ll]=0;
675 for(ll1=-1;ll1<boxs->nlist;ll1++) {
676 if(ll1==-1) {
677 mlo=sim->mols->topd;
678 mhi=sim->mols->nd;
679 mlist=sim->mols->dead; }
680 else {
681 mlo=0;
682 mhi=sim->mols->nl[ll1];
683 mlist=sim->mols->live[ll1]; }
684 for(m=mlo;m<mhi;m++) {
685 mptr=mlist[m];
686 if(mptr->ident>0) {
687 bptr=pos2box(sim,mptr->pos);
688 mptr->box=bptr;
689 ll=sim->mols->listlookup[mptr->ident][mptr->mstate];
690 bptr->nmol[ll]++; }}}
691
692 for(b=0;b<nbox;b++) {
693 bptr=blist[b];
694 for(ll=0;ll<boxs->nlist;ll++) {
695 mxml=bptr->nmol[ll];
696 bptr->nmol[ll]=0;
697 if(bptr->maxmol[ll]<mxml) {
698 er=expandbox(bptr,(int)(mxml*1.5-bptr->maxmol[ll]),ll);
699 if(er) return 1; }}}
700 for(ll1=0;ll1<boxs->nlist;ll1++) {
701 mlo=0;
702 mhi=sim->mols->nl[ll1];
703 mlist=sim->mols->live[ll1];
704 for(m=mlo;m<mhi;m++) {
705 mptr=mlist[m];
706 if(mptr->ident>0) {
707 ll=sim->mols->listlookup[mptr->ident][mptr->mstate];
708 bptr=mptr->box;
709 bptr->mol[ll][bptr->nmol[ll]++]=mptr; }}}}
710
711 return 0; }
712
713
714 /* boxesupdatelists */
boxesupdatelists(simptr sim)715 int boxesupdatelists(simptr sim) {
716 int dim,d,nbox,b,b2,w,er,nneigh,nwall;
717 int *side,*indx;
718 boxssptr boxs;
719 boxptr *blist,bptr;
720 double flt1,flt2,mpbox;
721 int *ntemp,*wptemp;
722
723 dim=sim->dim;
724 boxs=sim->boxs;
725
726 if(!boxs) { // create superstructure if needed
727 er=boxsetsize(sim,"molperbox",4);
728 if(er) return 1;
729 boxs=sim->boxs; }
730
731 if(sim->mols && sim->mols->condition<SCparams) return 2;
732 if(boxs->blist) { // box superstructure
733 boxesfree(boxs->blist,boxs->nbox,boxs->nlist);
734 boxs->nbox=0; }
735 side=boxs->side;
736 mpbox=boxs->mpbox;
737 if(mpbox<=0 && boxs->boxsize<=0) mpbox=5;
738 if(mpbox>0) {
739 flt1=systemvolume(sim);
740 flt2=(double)molcount(sim,-5,NULL,MSall,-1);
741 flt1=pow(flt2/mpbox/flt1,1.0/dim); }
742 else flt1=1.0/boxs->boxsize;
743 nbox=1;
744 for(d=0;d<dim;d++) {
745 boxs->min[d]=sim->wlist[2*d]->pos;
746 side[d]=(int)ceil((sim->wlist[2*d+1]->pos-sim->wlist[2*d]->pos)*flt1);
747 if(!side[d]) side[d]=1;
748 boxs->size[d]=(sim->wlist[2*d+1]->pos-sim->wlist[2*d]->pos)/side[d];
749 nbox*=side[d]; }
750 boxs->boxvol=1.0;
751 for(d=0;d<dim;d++) boxs->boxvol*=boxs->size[d];
752
753 boxs->nlist=sim->mols?sim->mols->nlist:0; // individual boxes
754 boxs->nbox=nbox;
755 blist=boxs->blist=boxesalloc(nbox,dim,boxs->nlist);
756 if(!blist) return 1;
757
758 for(b=0;b<nbox;b++) add2indxZV(b,blist[b]->indx,side,dim); // box->indx
759
760 if(sim->accur>=3) {
761 ntemp=allocZV(intpower(3,dim)-1); // neigh
762 wptemp=allocZV(intpower(3,dim)-1);
763 if(!ntemp || !wptemp) return 1;
764 for(b=0;b<nbox;b++) {
765 bptr=blist[b];
766 indx=bptr->indx;
767 if(sim->accur<6)
768 nneigh=neighborZV(indx,ntemp,side,dim,0,NULL,&bptr->midneigh);
769 else if(sim->accur<9) {
770 for(w=0;w<2*dim;w++) wptemp[w]=(sim->wlist[w]->type=='p' && sim->srfss==NULL);
771 nneigh=neighborZV(indx,ntemp,side,dim,6,wptemp,&bptr->midneigh); }
772 else {
773 for(w=0;w<2*dim;w++) wptemp[w]=(sim->wlist[w]->type=='p' && sim->srfss==NULL);
774 nneigh=neighborZV(indx,ntemp,side,dim,7,wptemp,&bptr->midneigh); }
775 if(nneigh==-1) return 1;
776 w=0;
777 if(sim->accur>=6) for(b2=0;b2<nneigh;b2++) w+=wptemp[b2];
778
779 if(nneigh && nneigh!=bptr->nneigh) {
780 free(bptr->neigh);
781 free(bptr->wpneigh);
782 bptr->wpneigh=NULL;
783 bptr->neigh=(boxptr*) calloc(nneigh,sizeof(boxptr));
784 if(!bptr->neigh) return 1;
785 if(w) {
786 bptr->wpneigh=allocZV(nneigh);
787 if(!bptr->wpneigh) return 1; }
788 bptr->nneigh=nneigh; }
789 if(nneigh)
790 for(b2=0;b2<nneigh;b2++) bptr->neigh[b2]=blist[ntemp[b2]];
791 if(w) {
792 for(b2=0;b2<nneigh;b2++) bptr->wpneigh[b2]=wptemp[b2]; }}
793
794 neighborZV(NULL,NULL,NULL,0,-1,NULL,NULL);
795 freeZV(ntemp);
796 freeZV(wptemp); }
797
798 for(b=0;b<nbox;b++) { // box->nwall, wlist
799 bptr=blist[b];
800 indx=bptr->indx;
801 nwall=0;
802 for(d=0;d<dim;d++) {
803 if(indx[d]==0) nwall++;
804 if(indx[d]==side[d]-1) nwall++; }
805 if(nwall && nwall!=bptr->nwall) {
806 free(bptr->wlist);
807 bptr->wlist=(wallptr*) calloc(nwall,sizeof(wallptr));
808 if(!bptr->wlist) return 1;
809 bptr->nwall=nwall; }
810 if(nwall) {
811 w=0;
812 for(d=0;d<dim;d++) {
813 if(indx[d]==0) bptr->wlist[w++]=sim->wlist[2*d];
814 if(indx[d]==side[d]-1) bptr->wlist[w++]=sim->wlist[2*d+1]; }}}
815
816 return 0; }
817
818
819 /* boxesupdate */
boxesupdate(simptr sim)820 int boxesupdate(simptr sim) {
821 int er;
822
823 if(sim->dim<=0 || !sim->wlist) return 3;
824
825 if(!sim->boxs || sim->boxs->condition<=SClists) {
826 er=boxesupdatelists(sim);
827 if(er) return er;
828 boxsetcondition(sim->boxs,SCparams,1); }
829
830 if(sim->boxs->condition==SCparams) {
831 er=boxesupdateparams(sim);
832 if(er) return er;
833 boxsetcondition(sim->boxs,SCok,1); }
834 return 0; }
835
836
837 /******************************************************************************/
838 /*************************** core simulation functions ************************/
839 /******************************************************************************/
840
841
842 /* line2nextbox */
line2nextbox(simptr sim,double * pt1,double * pt2,boxptr bptr)843 boxptr line2nextbox(simptr sim,double *pt1,double *pt2,boxptr bptr) {
844 int dim,d,d2,boxside,boxside2,adrs,z1[DIMMAX],*side,sum,flag;
845 double *size,*min,crsmin,edge,crs;
846
847 if(pos2box(sim,pt2)==bptr) return NULL;
848 dim=sim->dim;
849 size=sim->boxs->size;
850 side=sim->boxs->side;
851 min=sim->boxs->min;
852 crsmin=1.01;
853 boxside2=0;
854 d2=0;
855 flag=0;
856 for(d=0;d<dim;d++) {
857 z1[d]=bptr->indx[d];
858 if(pt2[d]!=pt1[d]) {
859 boxside=(pt2[d]>pt1[d])?1:0; // 1 for high side, 0 for low side
860 sum=bptr->indx[d]+boxside;
861 if(sum>0 && sum<side[d]) {
862 edge=min[d]+(double)sum*size[d]; // absolute location of potential edge crossing
863 crs=(edge-pt1[d])/(pt2[d]-pt1[d]); // relative position of potential edge crossing on line
864 if(crs<crsmin) {
865 crsmin=crs;
866 d2=d;
867 boxside2=boxside;
868 flag=0; }
869 else if(crs==crsmin) {
870 if(boxside==0 && boxside2==0 && (flag==0 || flag==2)) flag=2;
871 else flag=1; }}}}
872
873 if(flag) {
874 for(d=0;d<dim;d++) {
875 if(pt2[d]!=pt1[d]) {
876 boxside=(pt2[d]>pt1[d])?1:0;
877 sum=bptr->indx[d]+boxside;
878 if(sum>0 && sum<side[d]) {
879 edge=min[d]+(double)sum*size[d];
880 crs=(edge-pt1[d])/(pt2[d]-pt1[d]);
881 if(crs==crsmin && (boxside==1 || flag==2))
882 z1[d]+=boxside?1:-1; }}}
883 adrs=indx2addZV(z1,sim->boxs->side,dim);
884 return sim->boxs->blist[adrs]; }
885
886 if(crsmin==1.01) return NULL;
887 z1[d2]+=boxside2?1:-1;
888 adrs=indx2addZV(z1,sim->boxs->side,dim);
889 return sim->boxs->blist[adrs]; }
890
891
892 /* reassignmolecs */
reassignmolecs(simptr sim,int diffusing,int reborn)893 int reassignmolecs(simptr sim,int diffusing,int reborn) {
894 int m,nmol,m2,ll,b,s;
895 boxptr bptr1;
896 boxssptr boxss;
897 surfacessptr srfss;
898 moleculeptr mptr,*mlist,*mlist2;
899 surfaceptr srf;
900
901 if(!sim->mols) return 0;
902 boxss=sim->boxs;
903 srfss=sim->srfss;
904
905 if(!reborn) { // all molecules
906 for(ll=0;ll<sim->mols->nlist;ll++)
907 if(sim->mols->listtype[ll]==MLTsystem)
908 if(diffusing==0 || sim->mols->diffuselist[ll]==1) {
909 for(b=0;b<boxss->nbox;b++) // clear out box list
910 boxss->blist[b]->nmol[ll]=0;
911 if(srfss)
912 for(s=0;s<srfss->nsrf;s++) // clear out surface list
913 srfss->srflist[s]->nmol[ll]=0;
914 nmol=sim->mols->nl[ll];
915 mlist=sim->mols->live[ll];
916 for(m=0;m<nmol;m++) {
917 mptr=mlist[m];
918 bptr1=pos2box(sim,mptr->pos);
919 mptr->box=bptr1; // add to new box
920 if(bptr1->nmol[ll]==bptr1->maxmol[ll])
921 if(expandbox(bptr1,1+bptr1->nmol[ll],ll)) return 1;
922 bptr1->mol[ll][bptr1->nmol[ll]++]=mptr;
923 if(mptr->pnl) { // add to surface
924 srf=mptr->pnl->srf;
925 if(srf->nmol[ll]==srf->maxmol[ll])
926 if(surfexpandmollist(srf,2*srf->nmol[ll]+1,ll)) return 1;
927 srf->mol[ll][srf->nmol[ll]++]=mptr; }}}}
928
929 else { // reborn molecules only
930 for(ll=0;ll<sim->mols->nlist;ll++)
931 if(sim->mols->listtype[ll]==MLTsystem)
932 if(diffusing==0 || sim->mols->diffuselist[ll]==1) {
933 nmol=sim->mols->nl[ll];
934 mlist=sim->mols->live[ll];
935 for(m=sim->mols->topl[ll];m<nmol;m++) {
936 mptr=mlist[m];
937 bptr1=pos2box(sim,mptr->pos);
938 if(mptr->box!=bptr1) {
939 mlist2=mptr->box->mol[ll]; // remove from current box
940 for(m2=0;mlist2[m2]!=mptr;m2++);
941 mlist2[m2]=mlist2[--mptr->box->nmol[ll]];
942 mptr->box=bptr1; // add to new box
943 if(bptr1->nmol[ll]==bptr1->maxmol[ll])
944 if(expandbox(bptr1,1+bptr1->nmol[ll],ll)) return 1;
945 bptr1->mol[ll][bptr1->nmol[ll]++]=mptr; }
946 if(mptr->pnl) { // add to surface
947 srf=mptr->pnl->srf; // there is no check for prior listing on a surface because that is impossible
948 if(srf->nmol[ll]==srf->maxmol[ll])
949 if(surfexpandmollist(srf,2*srf->nmol[ll]+1,ll)) return 1;
950 srf->mol[ll][srf->nmol[ll]++]=mptr; }}}}
951 return 0; }
952
953
954
955