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 "random2.h"
11 #include "smoldyn.h"
12 #include "smoldynfuncs.h"
13
14 /******************************************************************************/
15 /************************************ Walls ***********************************/
16 /******************************************************************************/
17
18
19 /******************************************************************************/
20 /****************************** Local declarations ****************************/
21 /******************************************************************************/
22
23 // low level utilities
24
25 // memory management
26 wallptr wallalloc(void);
27 void wallfree(wallptr wptr);
28 wallptr *wallsalloc(int dim);
29
30 // data structure output
31
32 // structure setup
33
34 // core simulation functions
35
36
37 /******************************************************************************/
38 /****************************** low level utilities ***************************/
39 /******************************************************************************/
40
41 /* systemrandpos */
systemrandpos(simptr sim,double * pos)42 void systemrandpos(simptr sim,double *pos) {
43 int d;
44
45 for(d=0;d<sim->dim;d++) pos[d]=unirandOOD(sim->wlist[2*d]->pos,sim->wlist[2*d+1]->pos);
46 return; }
47
48
49 /* systemvolume */
systemvolume(simptr sim)50 double systemvolume(simptr sim) {
51 int d;
52 double vol;
53
54 vol=1.0;
55 for(d=0;d<sim->dim;d++) vol*=sim->wlist[2*d+1]->pos-sim->wlist[2*d]->pos;
56 return vol; }
57
58
59 /* systemcorners */
systemcorners(simptr sim,double * poslo,double * poshi)60 void systemcorners(simptr sim,double *poslo,double *poshi) {
61 int d;
62
63 for(d=0;d<sim->dim;d++) {
64 if(poslo) poslo[d]=sim->wlist[2*d]->pos;
65 if(poshi) poshi[d]=sim->wlist[2*d+1]->pos; }
66 return; }
67
68
69 /* systemcenter */
systemcenter(simptr sim,double * center)70 void systemcenter(simptr sim,double *center) {
71 int d;
72
73 for(d=0;d<sim->dim;d++)
74 center[d]=0.5*(sim->wlist[2*d]->pos+sim->wlist[2*d+1]->pos);
75 return; }
76
77
78 /* systemdiagonal */
systemdiagonal(simptr sim)79 double systemdiagonal(simptr sim) {
80 int d;
81 double diagonal;
82
83 diagonal=0;
84 for(d=0;d<sim->dim;d++)
85 diagonal+=(sim->wlist[2*d+1]->pos-sim->wlist[2*d]->pos)*(sim->wlist[2*d+1]->pos-sim->wlist[2*d]->pos);
86 return sqrt(diagonal); }
87
88
89 /* posinsystem */
posinsystem(simptr sim,double * pos)90 int posinsystem(simptr sim,double *pos) {
91 int d;
92
93 for(d=0;d<sim->dim;d++) {
94 if(pos[d]<sim->wlist[2*d]->pos) return 0;
95 if(pos[d]>sim->wlist[2*d+1]->pos) return 0; }
96 return 1; }
97
98
99 /* wallcalcdist2 */
wallcalcdist2(simptr sim,double * pos1,double * pos2,int wpcode,double * vect)100 double wallcalcdist2(simptr sim,double *pos1,double *pos2,int wpcode,double *vect) {
101 double dist2,syssize;
102 int d,dim;
103
104 dim=sim->dim;
105 dist2=0;
106 for(d=0;d<dim;d++) {
107 if((wpcode>>(2*d)&3)==0) // no wrapping
108 vect[d]=pos2[d]-pos1[d];
109 else if((wpcode>>(2*d)&3)==1) { // wrap towards low side, so pos1 is small and pos2 is big
110 syssize=sim->wlist[2*d+1]->pos-sim->wlist[2*d]->pos;
111 vect[d]=pos2[d]-pos1[d]-syssize; }
112 else { // wrap towards high side, so pos1 is big and pos2 is small
113 syssize=sim->wlist[2*d+1]->pos-sim->wlist[2*d]->pos;
114 vect[d]=pos2[d]-pos1[d]+syssize; }
115 dist2+=vect[d]*vect[d]; }
116 return dist2; }
117
118
119 /******************************************************************************/
120 /******************************* memory management ****************************/
121 /******************************************************************************/
122
123 /* wallalloc */
wallalloc(void)124 wallptr wallalloc(void) {
125 wallptr wptr;
126
127 CHECKMEM(wptr=(wallptr) malloc(sizeof(struct wallstruct)));
128 wptr->wdim=0;
129 wptr->side=0;
130 wptr->pos=0;
131 wptr->type='r';
132 wptr->opp=NULL;
133 return wptr;
134 failure:
135 simLog(NULL,10,"Unable to allocate memory in wallalloc");
136 return NULL; }
137
138
139 /* wallfree */
wallfree(wallptr wptr)140 void wallfree(wallptr wptr) {
141 if(!wptr) return;
142 free(wptr);
143 return; }
144
145
146 /* wallsalloc */
wallsalloc(int dim)147 wallptr *wallsalloc(int dim) {
148 int w,d;
149 wallptr *wlist;
150
151 if(dim<1) return NULL;
152 CHECKMEM(wlist=(wallptr *) calloc(2*dim,sizeof(wallptr)));
153 for(w=0;w<2*dim;w++) wlist[w]=NULL;
154 for(w=0;w<2*dim;w++)
155 CHECKMEM(wlist[w]=wallalloc());
156 for(d=0;d<dim;d++) {
157 wlist[2*d]->wdim=wlist[2*d+1]->wdim=d;
158 wlist[2*d]->side=0;
159 wlist[2*d+1]->side=1;
160 wlist[2*d]->pos=0;
161 wlist[2*d+1]->pos=1;
162 wlist[2*d]->type=wlist[2*d+1]->type='r';
163 wlist[2*d]->opp=wlist[2*d+1];
164 wlist[2*d+1]->opp=wlist[2*d]; }
165 return wlist;
166 failure:
167 wallsfree(wlist,dim);
168 simLog(NULL,10,"Unable to allocate memory in wallsalloc");
169 return NULL; }
170
171
172 /* wallsfree */
wallsfree(wallptr * wlist,int dim)173 void wallsfree(wallptr *wlist,int dim) {
174 if(!wlist||dim<1) return;
175 for(dim--;dim>=0;dim--) {
176 wallfree(wlist[2*dim+1]);
177 wallfree(wlist[2*dim]); }
178 free(wlist);
179 return; }
180
181
182 /******************************************************************************/
183 /***************************** data structure output **************************/
184 /******************************************************************************/
185
186 /* walloutput */
walloutput(simptr sim)187 void walloutput(simptr sim) {
188 int w,d,dim;
189 wallptr *wlist,wptr;
190 double vol,poslo[DIMMAX]={0},poshi[DIMMAX] = {0};
191 char dimchar;
192
193 dim=sim->dim;
194 wlist=sim->wlist;
195 simLog(sim,2,"WALL PARAMETERS\n");
196 if(!wlist) {
197 simLog(sim,2," No walls defined for simulation\n\n");
198 return; }
199
200 for(w=0;w<2*dim;w++) {
201 wptr=wlist[w];
202 if(wptr->wdim==0) dimchar='x';
203 else if(wptr->wdim==1) dimchar='y';
204 else if(wptr->wdim==2) dimchar='z';
205 else dimchar='?';
206 simLog(sim,2," wall %i: dimension %c, at %g",w,dimchar,wptr->pos);
207 if(sim->srfss) {
208 simLog(sim,1,", non-interacting because surfaces are defined");
209 simLog(sim,2,"\n"); }
210 else if(wptr->type=='r') simLog(sim,2,", reflecting\n");
211 else if(wptr->type=='p') simLog(sim,2,", periodic\n");
212 else if(wptr->type=='a') simLog(sim,2,", absorbing\n");
213 else if(wptr->type=='t') simLog(sim,2,", transparent\n");
214 if(wlist[w+1-2*(w%2)]!=wptr->opp) simLog(sim,10," ERROR: opposing wall is incorrect\n"); }
215
216 vol=systemvolume(sim);
217 if(dim==1) simLog(sim,2," system length: %g\n",vol);
218 else if(dim==2) simLog(sim,2," system area: %g\n",vol);
219 else simLog(sim,2," system volume: %g\n",vol);
220
221 systemcorners(sim,poslo,poshi);
222 simLog(sim,2," system corners: (%g",poslo[0]);
223 for(d=1;d<dim;d++) simLog(sim,2,",%g",poslo[d]);
224 simLog(sim,2,") and (%g",poshi[0]);
225 for(d=1;d<dim;d++) simLog(sim,2,",%g",poshi[d]);
226 simLog(sim,2,")\n");
227 simLog(sim,2,"\n");
228
229 return; }
230
231
232 /* writewalls */
writewalls(simptr sim,FILE * fptr)233 void writewalls(simptr sim,FILE *fptr) {
234 int d,dim;
235
236 fprintf(fptr,"# Boundary parameters\n");
237 dim=sim->dim;
238 for(d=0;d<dim;d++) {
239 fprintf(fptr,"low_wall %i %g %c\n",d,sim->wlist[2*d]->pos,sim->wlist[2*d]->type);
240 fprintf(fptr,"high_wall %i %g %c\n",d,sim->wlist[2*d+1]->pos,sim->wlist[2*d+1]->type); }
241 fprintf(fptr,"\n");
242 return; }
243
244
245 /* checkwallparams */
checkwallparams(simptr sim,int * warnptr)246 int checkwallparams(simptr sim,int *warnptr) {
247 int d,dim,warn,error;
248 wallptr *wlist;
249 double lowwall[DIMMAX],highwall[DIMMAX],syslen;
250
251 error=warn=0;
252 dim=sim->dim;
253 wlist=sim->wlist;
254
255 systemcorners(sim,lowwall,highwall);
256 syslen=0;
257 for(d=0;d<dim;d++) syslen+=(highwall[d]-lowwall[d])*(highwall[d]-lowwall[d]);
258 syslen=sqrt(syslen);
259 if(syslen<=0) {error++;simLog(sim,10," ERROR: Total system size is zero\n");}
260
261 for(d=0;d<dim;d++)
262 if(lowwall[d]>=highwall[d]) {
263 simLog(sim,10," ERROR: low_wall positions need to be smaller than high_wall positions");
264 error++; }
265
266 if(!sim->srfss) {
267 for(d=0;d<dim;d++) // check for asymmetric periodic boundary condition
268 if(wlist[2*d]->type=='p'&&wlist[2*d+1]->type!='p') {
269 simLog(sim,5," WARNING: only one wall on dimension %i has a periodic boundary condition\n",d);
270 warn++; }}
271
272 if(warnptr) *warnptr=warn;
273 return error; }
274
275
276 /******************************************************************************/
277 /******************************** structure setup *****************************/
278 /******************************************************************************/
279
280
281 /* walladd */
walladd(simptr sim,int d,int highside,double pos,char type)282 int walladd(simptr sim,int d,int highside,double pos,char type) {
283 if(!sim->wlist) {
284 if(!sim->dim) return 2;
285 sim->wlist=wallsalloc(sim->dim);
286 if(!sim->wlist) return 1; }
287 d=2*d+highside;
288 sim->wlist[d]->pos=pos;
289 sim->wlist[d]->type=type;
290 boxsetcondition(sim->boxs,SClists,0);
291 return 0; }
292
293
294 /* wallsettype */
wallsettype(simptr sim,int d,int highside,char type)295 int wallsettype(simptr sim,int d,int highside,char type) {
296 int d2;
297
298 if(!sim->wlist) return 1;
299 if(d<0)
300 for(d2=0;d2<sim->dim;d2++) {
301 if(highside<0) {
302 sim->wlist[2*d2]->type=type;
303 sim->wlist[2*d2+1]->type=type; }
304 else
305 sim->wlist[2*d2+highside]->type=type; }
306 else {
307 if(highside<0) {
308 sim->wlist[2*d]->type=type;
309 sim->wlist[2*d+1]->type=type; }
310 else
311 sim->wlist[2*d+highside]->type=type; }
312
313 boxsetcondition(sim->boxs,SClists,0);
314 return 0; }
315
316
317 /******************************************************************************/
318 /*************************** core simulation functions ************************/
319 /******************************************************************************/
320
321
322 /* checkwalls */
checkwalls(simptr sim,int ll,int reborn,boxptr bptr)323 int checkwalls(simptr sim,int ll,int reborn,boxptr bptr) {
324 int nmol,w,d,m;
325 moleculeptr *mlist;
326 double pos2,diff,difi,step,**difstep;
327 wallptr wptr;
328
329 if(sim->srfss) return 0;
330 if(bptr) {
331 nmol=bptr->nmol[ll];
332 mlist=bptr->mol[ll]; }
333 else {
334 nmol=sim->mols->nl[ll];
335 mlist=sim->mols->live[ll]; }
336 if(!reborn) m=0;
337 else if(reborn&&!bptr) m=sim->mols->topl[ll];
338 else {m=0;simLog(sim,10,"SMOLDYN ERROR: in checkwalls, both bptr and reborn are defined");}
339
340 for(w=0;w<2*sim->dim;w++) {
341 wptr=sim->wlist[w];
342 d=wptr->wdim;
343 if(wptr->type=='r'&&wptr->side==0) { // reflective
344 pos2=2*wptr->pos;
345 for(m=0;m<nmol;m++)
346 if(mlist[m]->pos[d]<wptr->pos) {
347 sim->eventcount[ETwall]++;
348 mlist[m]->pos[d]=pos2-mlist[m]->pos[d];}}
349 else if(wptr->type=='r') {
350 pos2=2*wptr->pos;
351 for(m=0;m<nmol;m++)
352 if(mlist[m]->pos[d]>wptr->pos) {
353 sim->eventcount[ETwall]++;
354 mlist[m]->pos[d]=pos2-mlist[m]->pos[d];}}
355 else if(wptr->type=='p'&&wptr->side==0) { // periodic
356 pos2=wptr->opp->pos-wptr->pos;
357 for(m=0;m<nmol;m++)
358 if(mlist[m]->pos[d]<wptr->pos) {
359 sim->eventcount[ETwall]++;
360 mlist[m]->pos[d]+=pos2;
361 mlist[m]->posoffset[d]-=pos2; }}
362 else if(wptr->type=='p') {
363 pos2=wptr->opp->pos-wptr->pos;
364 for(m=0;m<nmol;m++)
365 if(mlist[m]->pos[d]>wptr->pos) {
366 sim->eventcount[ETwall]++;
367 mlist[m]->pos[d]+=pos2;
368 mlist[m]->posoffset[d]-=pos2; }}
369 else if(wptr->type=='a') { // absorbing
370 difstep=sim->mols->difstep;
371 for(m=0;m<nmol;m++) {
372 diff=wptr->pos-mlist[m]->pos[d];
373 difi=wptr->pos-mlist[m]->posx[d];
374 step=difstep[mlist[m]->ident][MSsoln];
375 if((!(wptr->side)&&diff>0)||(wptr->side&&diff<0)||coinrandD(exp(-2*difi*diff/step/step))) {
376 sim->eventcount[ETwall]++;
377 molkill(sim,mlist[m],ll,-1); }}}}
378 sim->mols->touch++;
379 return 0; }
380
381
382
383