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