1 |
mmeineke |
821 |
#include <iostream> |
2 |
|
|
#include <stdio.h> |
3 |
|
|
#include <stdlib.h> |
4 |
|
|
#include <string.h> |
5 |
|
|
#include <math.h> |
6 |
|
|
|
7 |
|
|
#include "simError.h" |
8 |
|
|
#include "SimInfo.hpp" |
9 |
|
|
#include "ReadWrite.hpp" |
10 |
mmeineke |
825 |
#include "SimSetup.hpp" |
11 |
mmeineke |
821 |
|
12 |
|
|
#include "MoLocator.hpp" |
13 |
|
|
#include "latticeBuilder.hpp" |
14 |
|
|
|
15 |
mmeineke |
825 |
|
16 |
mmeineke |
821 |
#define VERSION_MAJOR 0 |
17 |
|
|
#define VERSION_MINOR 1 |
18 |
|
|
|
19 |
|
|
char *programName; /*the name of the program */ |
20 |
|
|
void usage(void); |
21 |
|
|
int buildLatticeBilayer( double aLat, |
22 |
|
|
double bLat, |
23 |
mmeineke |
825 |
double leafSpacing, |
24 |
|
|
char* waterName, |
25 |
|
|
char* lipidName); |
26 |
mmeineke |
821 |
using namespace std; |
27 |
mmeineke |
825 |
SimInfo* mainInfo; |
28 |
mmeineke |
821 |
|
29 |
|
|
int main(int argC,char* argV[]){ |
30 |
|
|
|
31 |
|
|
int i,j; // loop counters |
32 |
|
|
|
33 |
|
|
char* outPrefix; // the output prefix |
34 |
|
|
|
35 |
|
|
char* conversionCheck; |
36 |
|
|
bool conversionError; |
37 |
|
|
bool optionError; |
38 |
|
|
|
39 |
|
|
char currentFlag; // used in parsing the flags |
40 |
|
|
bool done = false; // multipurpose boolean |
41 |
|
|
bool havePrefix; // boolean for the output prefix |
42 |
|
|
|
43 |
|
|
char* lipidName; |
44 |
|
|
char* waterName; |
45 |
|
|
bool haveWaterName, haveLipidName; |
46 |
|
|
|
47 |
|
|
bool haveSpacing, haveAlat, haveBlat; |
48 |
|
|
double aLat, bLat, leafSpacing; |
49 |
|
|
|
50 |
|
|
char* inName; |
51 |
|
|
|
52 |
mmeineke |
825 |
SimSetup* simInit; |
53 |
|
|
|
54 |
mmeineke |
821 |
// first things first, all of the initializations |
55 |
|
|
|
56 |
|
|
fflush(stdout); |
57 |
|
|
srand48( 1337 ); // the random number generator. |
58 |
|
|
initSimError(); // the error handler |
59 |
|
|
|
60 |
|
|
outPrefix = NULL; |
61 |
|
|
inName = NULL; |
62 |
|
|
|
63 |
|
|
conversionError = false; |
64 |
|
|
optionError = false; |
65 |
|
|
|
66 |
|
|
havePrefix = false; |
67 |
|
|
haveAlat = false; |
68 |
|
|
haveBlat = false; |
69 |
|
|
haveSpacing = false; |
70 |
|
|
|
71 |
|
|
programName = argV[0]; /*save the program name in case we need it*/ |
72 |
|
|
|
73 |
|
|
for( i = 1; i < argC; i++){ |
74 |
|
|
|
75 |
|
|
if(argV[i][0] =='-'){ |
76 |
|
|
|
77 |
|
|
// parse the option |
78 |
|
|
|
79 |
|
|
if(argV[i][1] == '-' ){ |
80 |
|
|
|
81 |
|
|
// parse long word options |
82 |
|
|
|
83 |
|
|
if( !strcmp( argV[i], "--version") ){ |
84 |
|
|
|
85 |
|
|
printf("\n" |
86 |
|
|
"staticProps version %d.%d\n" |
87 |
|
|
"\n", |
88 |
|
|
VERSION_MAJOR, VERSION_MINOR ); |
89 |
|
|
exit(0); |
90 |
|
|
|
91 |
|
|
} |
92 |
|
|
|
93 |
|
|
else if( !strcmp( argV[i], "--help") ){ |
94 |
|
|
|
95 |
|
|
usage(); |
96 |
|
|
exit(0); |
97 |
|
|
} |
98 |
|
|
|
99 |
|
|
// anything else is an error |
100 |
|
|
|
101 |
|
|
else{ |
102 |
|
|
fprintf( stderr, |
103 |
|
|
"Invalid option \"%s\"\n", argV[i] ); |
104 |
|
|
usage(); |
105 |
|
|
exit(0); |
106 |
|
|
} |
107 |
|
|
} |
108 |
|
|
|
109 |
|
|
else{ |
110 |
|
|
|
111 |
|
|
// parse single character options |
112 |
|
|
|
113 |
|
|
done =0; |
114 |
|
|
j = 1; |
115 |
|
|
currentFlag = argV[i][j]; |
116 |
|
|
while( (currentFlag != '\0') && (!done) ){ |
117 |
|
|
|
118 |
|
|
switch(currentFlag){ |
119 |
|
|
|
120 |
|
|
case 'o': |
121 |
|
|
// -o <prefix> => the output prefix. |
122 |
|
|
|
123 |
|
|
j++; |
124 |
|
|
currentFlag = argV[i][j]; |
125 |
|
|
|
126 |
|
|
if( currentFlag != '\0' ) optionError = true; |
127 |
|
|
|
128 |
|
|
if( optionError ){ |
129 |
|
|
sprintf( painCave.errMsg, |
130 |
|
|
"\n" |
131 |
|
|
"The -o flag should end an option sequence.\n" |
132 |
|
|
" example: -r <outname> *NOT* -or <outname>\n" ); |
133 |
|
|
usage(); |
134 |
|
|
painCave.isFatal = 1; |
135 |
|
|
simError(); |
136 |
|
|
} |
137 |
|
|
|
138 |
|
|
i++; |
139 |
|
|
if( i>=argC ){ |
140 |
|
|
sprintf( painCave.errMsg, |
141 |
|
|
"\n" |
142 |
|
|
"not enough arguments for -o\n"); |
143 |
|
|
usage(); |
144 |
|
|
painCave.isFatal = 1; |
145 |
|
|
simError(); |
146 |
|
|
} |
147 |
|
|
|
148 |
|
|
outPrefix = argV[i]; |
149 |
|
|
if( outPrefix[0] == '-' ) optionError = true; |
150 |
|
|
|
151 |
|
|
if( optionError ){ |
152 |
|
|
sprintf( painCave.errMsg, |
153 |
|
|
"\n" |
154 |
|
|
"\"%s\" is not a valid out prefix/name.\n" |
155 |
|
|
"Out prefix/name should not begin with a dash.\n", |
156 |
|
|
outPrefix ); |
157 |
|
|
usage(); |
158 |
|
|
painCave.isFatal = 1; |
159 |
|
|
simError(); |
160 |
|
|
} |
161 |
|
|
|
162 |
|
|
havePrefix = true; |
163 |
|
|
done = true; |
164 |
|
|
break; |
165 |
|
|
|
166 |
|
|
case 'l': |
167 |
|
|
// -l <lipidName> => the lipid name. |
168 |
|
|
|
169 |
|
|
j++; |
170 |
|
|
currentFlag = argV[i][j]; |
171 |
|
|
|
172 |
|
|
if( currentFlag != '\0' ) optionError = true; |
173 |
|
|
|
174 |
|
|
if( optionError ){ |
175 |
|
|
sprintf( painCave.errMsg, |
176 |
|
|
"\n" |
177 |
|
|
"The -l flag should end an option sequence.\n" |
178 |
|
|
" example: -rl <lipidName> *NOT* -lr <lipidName>\n" ); |
179 |
|
|
usage(); |
180 |
|
|
painCave.isFatal = 1; |
181 |
|
|
simError(); |
182 |
|
|
} |
183 |
|
|
|
184 |
|
|
i++; |
185 |
|
|
if( i>=argC ){ |
186 |
|
|
sprintf( painCave.errMsg, |
187 |
|
|
"\n" |
188 |
|
|
"not enough arguments for -l\n"); |
189 |
|
|
usage(); |
190 |
|
|
painCave.isFatal = 1; |
191 |
|
|
simError(); |
192 |
|
|
} |
193 |
|
|
|
194 |
|
|
lipidName = argV[i]; |
195 |
|
|
if( lipidName[0] == '-' ) optionError = true; |
196 |
|
|
|
197 |
|
|
if( optionError ){ |
198 |
|
|
sprintf( painCave.errMsg, |
199 |
|
|
"\n" |
200 |
|
|
"\"%s\" is not a valid lipidName.\n" |
201 |
|
|
"lipidName should not begin with a dash.\n", |
202 |
|
|
lipidName ); |
203 |
|
|
usage(); |
204 |
|
|
painCave.isFatal = 1; |
205 |
|
|
simError(); |
206 |
|
|
} |
207 |
|
|
|
208 |
|
|
haveLipidName = true; |
209 |
|
|
done = true; |
210 |
|
|
break; |
211 |
|
|
|
212 |
|
|
case 'w': |
213 |
|
|
// -w <waterName> => the water name. |
214 |
|
|
|
215 |
|
|
j++; |
216 |
|
|
currentFlag = argV[i][j]; |
217 |
|
|
|
218 |
|
|
if( currentFlag != '\0' ) optionError = true; |
219 |
|
|
|
220 |
|
|
if( optionError ){ |
221 |
|
|
sprintf( painCave.errMsg, |
222 |
|
|
"\n" |
223 |
|
|
"The -w flag should end an option sequence.\n" |
224 |
|
|
" example: -rw <waterName> *NOT* -lw <waterName>\n" ); |
225 |
|
|
usage(); |
226 |
|
|
painCave.isFatal = 1; |
227 |
|
|
simError(); |
228 |
|
|
} |
229 |
|
|
|
230 |
|
|
i++; |
231 |
|
|
if( i>=argC ){ |
232 |
|
|
sprintf( painCave.errMsg, |
233 |
|
|
"\n" |
234 |
|
|
"not enough arguments for -w\n"); |
235 |
|
|
usage(); |
236 |
|
|
painCave.isFatal = 1; |
237 |
|
|
simError(); |
238 |
|
|
} |
239 |
|
|
|
240 |
|
|
waterName = argV[i]; |
241 |
|
|
if( waterName[0] == '-' ) optionError = true; |
242 |
|
|
|
243 |
|
|
if( optionError ){ |
244 |
|
|
sprintf( painCave.errMsg, |
245 |
|
|
"\n" |
246 |
|
|
"\"%s\" is not a valid waterName.\n" |
247 |
|
|
"waterName should not begin with a dash.\n", |
248 |
|
|
waterName ); |
249 |
|
|
usage(); |
250 |
|
|
painCave.isFatal = 1; |
251 |
|
|
simError(); |
252 |
|
|
} |
253 |
|
|
|
254 |
|
|
haveWaterName = true; |
255 |
|
|
done = true; |
256 |
|
|
break; |
257 |
|
|
|
258 |
|
|
|
259 |
|
|
case 'h': |
260 |
|
|
// -h <double> set <double> to the hex lattice spacing |
261 |
|
|
|
262 |
|
|
haveAlat = true; |
263 |
|
|
haveBlat = true; |
264 |
|
|
j++; |
265 |
|
|
currentFlag = argV[i][j]; |
266 |
|
|
|
267 |
|
|
if( currentFlag != '\0' ) optionError = true; |
268 |
|
|
|
269 |
|
|
if( optionError ){ |
270 |
|
|
sprintf( painCave.errMsg, |
271 |
|
|
"\n" |
272 |
|
|
"The -h flag should end an option sequence.\n" |
273 |
|
|
" example: -sh <double> *NOT* -hs <double>\n" ); |
274 |
|
|
usage(); |
275 |
|
|
painCave.isFatal = 1; |
276 |
|
|
simError(); |
277 |
|
|
} |
278 |
|
|
|
279 |
|
|
i++; |
280 |
|
|
if( i>=argC ){ |
281 |
|
|
sprintf( painCave.errMsg, |
282 |
|
|
"\n" |
283 |
|
|
"not enough arguments for -h\n"); |
284 |
|
|
usage(); |
285 |
|
|
painCave.isFatal = 1; |
286 |
|
|
simError(); |
287 |
|
|
} |
288 |
|
|
|
289 |
mmeineke |
825 |
|
290 |
|
|
bLat = atof( argV[i] ); |
291 |
|
|
// bLat = strtod( argV[i], &conversionCheck); |
292 |
|
|
// if( conversionCheck == argV[i] ) conversionError = true; |
293 |
|
|
// if( *conversionCheck != '\0' ) conversionError = true; |
294 |
mmeineke |
821 |
|
295 |
|
|
if( conversionError ){ |
296 |
|
|
sprintf( painCave.errMsg, |
297 |
|
|
"Error converting \"%s\" to a double for \"h\".\n", |
298 |
|
|
argV[i] ); |
299 |
|
|
usage(); |
300 |
|
|
painCave.isFatal = 1; |
301 |
|
|
simError(); |
302 |
|
|
} |
303 |
|
|
|
304 |
|
|
aLat = sqrt( 3.0 ) * bLat; |
305 |
|
|
|
306 |
|
|
done = true; |
307 |
|
|
|
308 |
|
|
break; |
309 |
|
|
|
310 |
|
|
case 'a': |
311 |
|
|
// -a <double> set <double> to the aLat |
312 |
|
|
|
313 |
|
|
haveAlat = true; |
314 |
|
|
j++; |
315 |
|
|
currentFlag = argV[i][j]; |
316 |
|
|
|
317 |
|
|
if( currentFlag != '\0' ) optionError = true; |
318 |
|
|
|
319 |
|
|
if( optionError ){ |
320 |
|
|
sprintf( painCave.errMsg, |
321 |
|
|
"\n" |
322 |
|
|
"The -a flag should end an option sequence.\n" |
323 |
|
|
" example: -sa <double> *NOT* -as <double>\n" ); |
324 |
|
|
usage(); |
325 |
|
|
painCave.isFatal = 1; |
326 |
|
|
simError(); |
327 |
|
|
} |
328 |
|
|
|
329 |
|
|
i++; |
330 |
|
|
if( i>=argC ){ |
331 |
|
|
sprintf( painCave.errMsg, |
332 |
|
|
"\n" |
333 |
|
|
"not enough arguments for -a\n"); |
334 |
|
|
usage(); |
335 |
|
|
painCave.isFatal = 1; |
336 |
|
|
simError(); |
337 |
|
|
} |
338 |
|
|
|
339 |
mmeineke |
825 |
aLat = atof( argV[i] ); |
340 |
|
|
// aLat = strtod( argV[i], &conversionCheck); |
341 |
|
|
// if( conversionCheck == argV[i] ) conversionError = true; |
342 |
|
|
// if( *conversionCheck != '\0' ) conversionError = true; |
343 |
mmeineke |
821 |
|
344 |
|
|
if( conversionError ){ |
345 |
|
|
sprintf( painCave.errMsg, |
346 |
|
|
"Error converting \"%s\" to a double for \"a\".\n", |
347 |
|
|
argV[i] ); |
348 |
|
|
usage(); |
349 |
|
|
painCave.isFatal = 1; |
350 |
|
|
simError(); |
351 |
|
|
} |
352 |
|
|
|
353 |
|
|
done = true; |
354 |
|
|
|
355 |
|
|
break; |
356 |
|
|
|
357 |
|
|
case 'b': |
358 |
|
|
// -b <double> set <double> to the bLat |
359 |
|
|
|
360 |
|
|
haveBlat = true; |
361 |
|
|
j++; |
362 |
|
|
currentFlag = argV[i][j]; |
363 |
|
|
|
364 |
|
|
if( currentFlag != '\0' ) optionError = true; |
365 |
|
|
|
366 |
|
|
if( optionError ){ |
367 |
|
|
sprintf( painCave.errMsg, |
368 |
|
|
"\n" |
369 |
|
|
"The -b flag should end an option sequence.\n" |
370 |
|
|
" example: -sb <double> *NOT* -bs <double>\n" ); |
371 |
|
|
usage(); |
372 |
|
|
painCave.isFatal = 1; |
373 |
|
|
simError(); |
374 |
|
|
} |
375 |
|
|
|
376 |
|
|
i++; |
377 |
|
|
if( i>=argC ){ |
378 |
|
|
sprintf( painCave.errMsg, |
379 |
|
|
"\n" |
380 |
|
|
"not enough arguments for -b\n"); |
381 |
|
|
usage(); |
382 |
|
|
painCave.isFatal = 1; |
383 |
|
|
simError(); |
384 |
|
|
} |
385 |
|
|
|
386 |
mmeineke |
825 |
bLat = atof( argV[i] ); |
387 |
|
|
|
388 |
|
|
// bLat = strtod( argV[i], &conversionCheck); |
389 |
|
|
// if( conversionCheck == argV[i] ) conversionError = true; |
390 |
|
|
// if( *conversionCheck != '\0' ) conversionError = true; |
391 |
mmeineke |
821 |
|
392 |
|
|
if( conversionError ){ |
393 |
|
|
sprintf( painCave.errMsg, |
394 |
|
|
"Error converting \"%s\" to a double for \"a\".\n", |
395 |
|
|
argV[i] ); |
396 |
|
|
usage(); |
397 |
|
|
painCave.isFatal = 1; |
398 |
|
|
simError(); |
399 |
|
|
} |
400 |
|
|
|
401 |
|
|
done = true; |
402 |
|
|
|
403 |
|
|
break; |
404 |
|
|
|
405 |
|
|
case 's': |
406 |
|
|
// -s <double> set <double> to the leafSpacing |
407 |
|
|
|
408 |
|
|
haveSpacing = true; |
409 |
|
|
j++; |
410 |
|
|
currentFlag = argV[i][j]; |
411 |
|
|
|
412 |
|
|
if( currentFlag != '\0' ) optionError = true; |
413 |
|
|
|
414 |
|
|
if( optionError ){ |
415 |
|
|
sprintf( painCave.errMsg, |
416 |
|
|
"\n" |
417 |
|
|
"The -s flag should end an option sequence.\n" |
418 |
|
|
" example: -rs <double> *NOT* -sr <double>\n" ); |
419 |
|
|
usage(); |
420 |
|
|
painCave.isFatal = 1; |
421 |
|
|
simError(); |
422 |
|
|
} |
423 |
|
|
|
424 |
|
|
i++; |
425 |
|
|
if( i>=argC ){ |
426 |
|
|
sprintf( painCave.errMsg, |
427 |
|
|
"\n" |
428 |
|
|
"not enough arguments for -s\n"); |
429 |
|
|
usage(); |
430 |
|
|
painCave.isFatal = 1; |
431 |
|
|
simError(); |
432 |
|
|
} |
433 |
|
|
|
434 |
mmeineke |
825 |
leafSpacing = atof( argV[i] ); |
435 |
|
|
|
436 |
|
|
// leafSpacing = strtod( argV[i], &conversionCheck); |
437 |
|
|
// if( conversionCheck == argV[i] ) conversionError = true; |
438 |
|
|
// if( *conversionCheck != '\0' ) conversionError = true; |
439 |
mmeineke |
821 |
|
440 |
|
|
if( conversionError ){ |
441 |
|
|
sprintf( painCave.errMsg, |
442 |
|
|
"Error converting \"%s\" to a double for \"s\".\n", |
443 |
|
|
argV[i] ); |
444 |
|
|
usage(); |
445 |
|
|
painCave.isFatal = 1; |
446 |
|
|
simError(); |
447 |
|
|
} |
448 |
|
|
|
449 |
|
|
done = true; |
450 |
|
|
|
451 |
|
|
break; |
452 |
|
|
|
453 |
|
|
default: |
454 |
|
|
|
455 |
|
|
sprintf(painCave.errMsg, |
456 |
|
|
"\n" |
457 |
|
|
"Bad option \"-%c\"\n", currentFlag); |
458 |
|
|
usage(); |
459 |
|
|
painCave.isFatal = 1; |
460 |
|
|
simError(); |
461 |
|
|
} |
462 |
|
|
j++; |
463 |
|
|
currentFlag = argV[i][j]; |
464 |
|
|
} |
465 |
|
|
} |
466 |
|
|
} |
467 |
|
|
|
468 |
|
|
else{ |
469 |
|
|
|
470 |
|
|
if( inName != NULL ){ |
471 |
|
|
sprintf( painCave.errMsg, |
472 |
|
|
"Error at \"%s\", program does not currently support\n" |
473 |
|
|
"more than one input bass file.\n" |
474 |
|
|
"\n", |
475 |
|
|
argV[i]); |
476 |
|
|
usage(); |
477 |
|
|
painCave.isFatal = 1; |
478 |
|
|
simError(); |
479 |
|
|
} |
480 |
|
|
|
481 |
|
|
inName = argV[i]; |
482 |
|
|
|
483 |
|
|
} |
484 |
|
|
} |
485 |
|
|
|
486 |
|
|
if( inName == NULL ){ |
487 |
|
|
sprintf( painCave.errMsg, |
488 |
|
|
"Error, bass file is needed to run.\n" ); |
489 |
|
|
usage(); |
490 |
|
|
painCave.isFatal = 1; |
491 |
|
|
simError(); |
492 |
|
|
} |
493 |
|
|
|
494 |
|
|
// if no output prefix is given default to "donkey". |
495 |
|
|
|
496 |
|
|
if( !havePrefix ){ |
497 |
|
|
outPrefix = strdup( "donkey" ); |
498 |
|
|
} |
499 |
|
|
|
500 |
|
|
|
501 |
|
|
if( !haveWaterName ){ |
502 |
|
|
sprintf( painCave.errMsg, |
503 |
|
|
"Error, the water name is needed to run.\n" |
504 |
|
|
); |
505 |
|
|
usage(); |
506 |
|
|
painCave.isFatal = 1; |
507 |
|
|
simError(); |
508 |
|
|
} |
509 |
|
|
|
510 |
|
|
if( !haveLipidName ){ |
511 |
|
|
sprintf( painCave.errMsg, |
512 |
|
|
"Error, the lipid name is needed to run.\n" |
513 |
|
|
); |
514 |
|
|
usage(); |
515 |
|
|
painCave.isFatal = 1; |
516 |
|
|
simError(); |
517 |
|
|
} |
518 |
|
|
|
519 |
|
|
if( !haveAlat ){ |
520 |
|
|
sprintf( painCave.errMsg, |
521 |
|
|
"Error, the hexagonal lattice parameter, a, is needed to run.\n" |
522 |
|
|
); |
523 |
|
|
usage(); |
524 |
|
|
painCave.isFatal = 1; |
525 |
|
|
simError(); |
526 |
|
|
} |
527 |
|
|
|
528 |
|
|
if( !haveBlat ){ |
529 |
|
|
sprintf( painCave.errMsg, |
530 |
|
|
"Error, the hexagonal lattice parameter, b, is needed to run.\n" |
531 |
|
|
); |
532 |
|
|
usage(); |
533 |
|
|
painCave.isFatal = 1; |
534 |
|
|
simError(); |
535 |
|
|
} |
536 |
|
|
|
537 |
|
|
|
538 |
|
|
if( !haveSpacing ){ |
539 |
|
|
sprintf( painCave.errMsg, |
540 |
|
|
"Error, the leaf spaceing is needed to run.\n" |
541 |
|
|
); |
542 |
|
|
usage(); |
543 |
|
|
painCave.isFatal = 1; |
544 |
|
|
simError(); |
545 |
|
|
} |
546 |
|
|
|
547 |
mmeineke |
825 |
mainInfo = new SimInfo(); |
548 |
|
|
simInit = new SimSetup(); |
549 |
|
|
simInit->setSimInfo( mainInfo ); |
550 |
|
|
simInit->suspendInit(); |
551 |
|
|
simInit->parseFile( inName ); |
552 |
mmeineke |
874 |
simInit->createSim(); |
553 |
mmeineke |
821 |
|
554 |
mmeineke |
825 |
delete simInit; |
555 |
mmeineke |
821 |
|
556 |
mmeineke |
825 |
sprintf( mainInfo->statusName, "%s.stat", outPrefix ); |
557 |
|
|
sprintf( mainInfo->sampleName, "%s.dump", outPrefix ); |
558 |
|
|
sprintf( mainInfo->finalName, "%s.init", outPrefix ); |
559 |
mmeineke |
821 |
|
560 |
mmeineke |
825 |
buildLatticeBilayer( aLat, bLat, leafSpacing, waterName, lipidName ); |
561 |
|
|
|
562 |
mmeineke |
821 |
return 0; |
563 |
|
|
} |
564 |
|
|
|
565 |
|
|
int buildLatticeBilayer(double aLat, |
566 |
|
|
double bLat, |
567 |
mmeineke |
825 |
double leafSpacing, |
568 |
|
|
char* waterName, |
569 |
|
|
char* lipidName){ |
570 |
mmeineke |
821 |
|
571 |
|
|
typedef struct{ |
572 |
|
|
double rot[3][3]; |
573 |
|
|
double pos[3]; |
574 |
|
|
} coord; |
575 |
|
|
|
576 |
|
|
const double waterRho = 0.0334; // number density per cubic angstrom |
577 |
|
|
const double waterVol = 4.0 / waterRho; // volume occupied by 4 waters |
578 |
|
|
|
579 |
|
|
double waterCell[3]; |
580 |
|
|
|
581 |
|
|
double *posX, *posY, *posZ; |
582 |
|
|
double pos[3], posA[3], posB[3]; |
583 |
|
|
|
584 |
mmeineke |
874 |
const double waterFudge = 6.0; |
585 |
mmeineke |
821 |
|
586 |
|
|
int i,j,k,l; |
587 |
|
|
int nAtoms, atomIndex, molIndex, molID; |
588 |
|
|
int* molSeq; |
589 |
|
|
int* molMap; |
590 |
|
|
int* molStart; |
591 |
|
|
int testTot, done; |
592 |
|
|
int nCells, nCellsX, nCellsY, nCellsZ; |
593 |
|
|
int nx, ny; |
594 |
|
|
double boxX, boxY, boxZ; |
595 |
mmeineke |
825 |
double theta, phi, psi; |
596 |
mmeineke |
821 |
int which; |
597 |
|
|
int targetWaters; |
598 |
|
|
|
599 |
|
|
Atom** atoms; |
600 |
|
|
SimInfo* testInfo; |
601 |
|
|
coord testSite; |
602 |
|
|
SimState* theConfig; |
603 |
|
|
DumpWriter* writer; |
604 |
|
|
|
605 |
|
|
MoleculeStamp* lipidStamp; |
606 |
|
|
MoleculeStamp* waterStamp; |
607 |
|
|
MoLocator *lipidLocate; |
608 |
|
|
MoLocator *waterLocate; |
609 |
|
|
int foundLipid, foundWater; |
610 |
|
|
int nLipids, lipidNatoms, nWaters, waterNatoms; |
611 |
|
|
int targetNlipids, targetNwaters; |
612 |
|
|
double targetWaterLipidRatio; |
613 |
|
|
double maxZ, minZ, zHeight; |
614 |
|
|
|
615 |
mmeineke |
825 |
molStart = NULL; |
616 |
|
|
|
617 |
mmeineke |
821 |
// set the the lipidStamp |
618 |
|
|
|
619 |
|
|
foundLipid = 0; |
620 |
|
|
foundWater = 0; |
621 |
mmeineke |
825 |
|
622 |
|
|
for(i=0; i<mainInfo->nComponents; i++){ |
623 |
|
|
|
624 |
|
|
if( !strcmp( mainInfo->compStamps[i]->getID(), lipidName ) ){ |
625 |
mmeineke |
821 |
|
626 |
|
|
foundLipid = 1; |
627 |
mmeineke |
825 |
lipidStamp = mainInfo->compStamps[i]; |
628 |
|
|
targetNlipids = mainInfo->componentsNmol[i]; |
629 |
mmeineke |
821 |
lipidNatoms = lipidStamp->getNAtoms(); |
630 |
|
|
} |
631 |
mmeineke |
825 |
if( !strcmp( mainInfo->compStamps[i]->getID(), waterName ) ){ |
632 |
mmeineke |
821 |
|
633 |
|
|
foundWater = 1; |
634 |
|
|
|
635 |
mmeineke |
825 |
waterStamp = mainInfo->compStamps[i]; |
636 |
|
|
targetNwaters = mainInfo->componentsNmol[i]; |
637 |
mmeineke |
821 |
waterNatoms = waterStamp->getNAtoms(); |
638 |
|
|
} |
639 |
|
|
} |
640 |
|
|
if( !foundLipid ){ |
641 |
|
|
sprintf(painCave.errMsg, |
642 |
mmeineke |
825 |
"latticeBilayer error: Could not find lipid \"%s\" in the bass file.\n", |
643 |
|
|
lipidName ); |
644 |
mmeineke |
821 |
painCave.isFatal = 1; |
645 |
|
|
simError(); |
646 |
|
|
} |
647 |
|
|
if( !foundWater ){ |
648 |
|
|
sprintf(painCave.errMsg, |
649 |
mmeineke |
825 |
"latticeBilayer error: Could not find solvent \"%s\" in the bass file.\n", |
650 |
|
|
waterName ); |
651 |
mmeineke |
821 |
painCave.isFatal = 1; |
652 |
|
|
simError(); |
653 |
|
|
} |
654 |
|
|
|
655 |
|
|
//create the Molocator arrays |
656 |
|
|
|
657 |
|
|
lipidLocate = new MoLocator( lipidStamp ); |
658 |
|
|
waterLocate = new MoLocator( waterStamp ); |
659 |
|
|
|
660 |
|
|
// gather info about the lipid |
661 |
|
|
|
662 |
|
|
testInfo = new SimInfo(); |
663 |
|
|
testInfo->n_atoms = lipidNatoms; |
664 |
|
|
theConfig = testInfo->getConfiguration(); |
665 |
|
|
theConfig->createArrays( lipidNatoms ); |
666 |
|
|
testInfo->atoms = new Atom*[lipidNatoms]; |
667 |
|
|
atoms = testInfo->atoms; |
668 |
|
|
|
669 |
|
|
testSite.pos[0] = 0.0; |
670 |
|
|
testSite.pos[1] = 0.0; |
671 |
|
|
testSite.pos[2] = 0.0; |
672 |
|
|
|
673 |
mmeineke |
825 |
theta = 0.0; |
674 |
|
|
phi = 0.0; |
675 |
|
|
psi = 0.0; |
676 |
|
|
|
677 |
|
|
getEulerRot(theta, phi, psi, testSite.rot ); |
678 |
mmeineke |
821 |
lipidLocate->placeMol( testSite.pos, testSite.rot, atoms, 0, theConfig ); |
679 |
|
|
|
680 |
|
|
minZ = 0.0; |
681 |
|
|
maxZ = 0.0; |
682 |
|
|
double myPos[3]; |
683 |
|
|
for(i=0;i<lipidNatoms;i++){ |
684 |
|
|
atoms[i]->getPos( myPos ); |
685 |
|
|
minZ = (minZ > myPos[2]) ? myPos[2] : minZ; |
686 |
|
|
maxZ = (maxZ < myPos[2]) ? myPos[2] : maxZ; |
687 |
|
|
} |
688 |
|
|
zHeight = maxZ - minZ; |
689 |
|
|
|
690 |
mmeineke |
874 |
std::cerr << "aLat = " << aLat << "; bLat = " << bLat << "\n"; |
691 |
mmeineke |
821 |
|
692 |
mmeineke |
874 |
nCells = (int)ceil( sqrt( (double)targetNlipids * bLat / (4.0 * aLat) )); |
693 |
|
|
|
694 |
mmeineke |
821 |
nx = nCells; |
695 |
|
|
ny = (int) ((double)nCells * aLat / bLat); |
696 |
|
|
|
697 |
mmeineke |
874 |
std::cerr << "nx = " << nx << "; ny = " << ny << "\n"; |
698 |
|
|
|
699 |
mmeineke |
821 |
boxX = nx * aLat; |
700 |
|
|
boxY = ny * bLat; |
701 |
|
|
|
702 |
|
|
nLipids = 4 * nx * ny; |
703 |
|
|
coord* lipidSites = new coord[nLipids]; |
704 |
|
|
|
705 |
mmeineke |
825 |
phi = 0.0; |
706 |
|
|
psi = 0.0; |
707 |
mmeineke |
821 |
|
708 |
|
|
which = 0; |
709 |
|
|
|
710 |
|
|
for (i = 0; i < nx; i++) { |
711 |
|
|
for (j = 0; j < ny; j++ ) { |
712 |
|
|
for (k = 0; k < 2; k++) { |
713 |
|
|
|
714 |
|
|
lipidSites[which].pos[0] = (double)i * aLat; |
715 |
|
|
lipidSites[which].pos[1] = (double)j * bLat; |
716 |
mmeineke |
825 |
lipidSites[which].pos[2] = (2.0* (double)k - 1.0) * |
717 |
|
|
((leafSpacing / 2.0) - maxZ); |
718 |
mmeineke |
821 |
|
719 |
mmeineke |
825 |
theta = (1.0 - (double)k) * M_PI; |
720 |
mmeineke |
821 |
|
721 |
mmeineke |
825 |
getEulerRot( theta, phi, psi, lipidSites[which].rot ); |
722 |
mmeineke |
821 |
|
723 |
|
|
which++; |
724 |
|
|
|
725 |
|
|
lipidSites[which].pos[0] = aLat * ((double)i + 0.5); |
726 |
|
|
lipidSites[which].pos[1] = bLat * ((double)j + 0.5); |
727 |
mmeineke |
825 |
lipidSites[which].pos[2] = (2.0* (double)k - 1.0) * |
728 |
|
|
((leafSpacing / 2.0) - maxZ); |
729 |
mmeineke |
821 |
|
730 |
mmeineke |
825 |
theta = (1.0 - (double)k) * M_PI; |
731 |
mmeineke |
821 |
|
732 |
mmeineke |
825 |
getEulerRot( theta, phi, psi, lipidSites[which].rot ); |
733 |
mmeineke |
821 |
|
734 |
|
|
which++; |
735 |
|
|
} |
736 |
|
|
} |
737 |
|
|
} |
738 |
|
|
|
739 |
|
|
targetWaterLipidRatio = (double)targetNwaters / (double)targetNlipids; |
740 |
|
|
|
741 |
|
|
targetWaters = targetWaterLipidRatio * nLipids; |
742 |
|
|
|
743 |
|
|
// guess the size of the water box |
744 |
|
|
|
745 |
|
|
nCellsX = (int)ceil(boxX / pow(waterVol, ( 1.0 / 3.0 )) ); |
746 |
|
|
nCellsY = (int)ceil(boxY / pow(waterVol, ( 1.0 / 3.0 )) ); |
747 |
|
|
|
748 |
|
|
done = 0; |
749 |
|
|
nCellsZ = 0; |
750 |
|
|
while( !done ){ |
751 |
|
|
|
752 |
|
|
nCellsZ++; |
753 |
|
|
testTot = 4 * nCellsX * nCellsY * nCellsZ; |
754 |
|
|
|
755 |
|
|
if( testTot >= targetWaters ) done = 1; |
756 |
|
|
} |
757 |
|
|
|
758 |
|
|
nWaters = nCellsX * nCellsY * nCellsZ * 4; |
759 |
|
|
|
760 |
mmeineke |
825 |
// calc current system size; |
761 |
|
|
|
762 |
|
|
nAtoms = 0; |
763 |
|
|
molIndex = 0; |
764 |
|
|
if(molStart != NULL ) delete[] molStart; |
765 |
|
|
molStart = new int[nLipids]; |
766 |
|
|
|
767 |
|
|
for(j=0; j<nLipids; j++){ |
768 |
|
|
molStart[molIndex] = nAtoms; |
769 |
|
|
molIndex++; |
770 |
|
|
nAtoms += lipidNatoms; |
771 |
|
|
} |
772 |
|
|
|
773 |
|
|
testInfo->n_atoms = nAtoms; |
774 |
|
|
theConfig = testInfo->getConfiguration(); |
775 |
|
|
theConfig->destroyArrays(); |
776 |
|
|
theConfig->createArrays( nAtoms ); |
777 |
|
|
testInfo->atoms = new Atom*[nAtoms]; |
778 |
|
|
atoms = testInfo->atoms; |
779 |
|
|
|
780 |
|
|
molIndex = 0; |
781 |
|
|
for(i=0; i<nLipids; i++ ){ |
782 |
|
|
lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms, |
783 |
|
|
molStart[molIndex], theConfig ); |
784 |
|
|
molIndex++; |
785 |
|
|
} |
786 |
|
|
|
787 |
|
|
atoms[0]->getPos( myPos ); |
788 |
|
|
|
789 |
|
|
maxZ = myPos[2]; |
790 |
|
|
minZ = myPos[2]; |
791 |
|
|
|
792 |
|
|
for(i=0;i<nAtoms;i++){ |
793 |
|
|
atoms[i]->getPos( myPos ); |
794 |
|
|
|
795 |
|
|
minZ = (minZ > myPos[2]) ? myPos[2] : minZ; |
796 |
|
|
maxZ = (maxZ < myPos[2]) ? myPos[2] : maxZ; |
797 |
|
|
} |
798 |
|
|
|
799 |
|
|
boxZ = (maxZ - minZ)+2.0; |
800 |
|
|
|
801 |
|
|
double centerX, centerY, centerZ; |
802 |
|
|
|
803 |
mmeineke |
943 |
centerX = (boxX / 2.0); |
804 |
|
|
centerY = (boxY / 2.0); |
805 |
mmeineke |
825 |
centerZ = ((maxZ - minZ) / 2.0) + minZ; |
806 |
|
|
|
807 |
|
|
// set up water coordinates |
808 |
|
|
|
809 |
mmeineke |
821 |
coord* waterSites = new coord[nWaters]; |
810 |
|
|
|
811 |
|
|
waterCell[0] = boxX / nCellsX; |
812 |
|
|
waterCell[1] = boxY / nCellsY; |
813 |
|
|
waterCell[2] = 4.0 / (waterRho * waterCell[0] * waterCell[1]); |
814 |
|
|
|
815 |
|
|
Lattice *myORTHO; |
816 |
|
|
myORTHO = new Lattice( ORTHORHOMBIC_LATTICE_TYPE, waterCell); |
817 |
mmeineke |
825 |
myORTHO->setStartZ( maxZ + waterFudge); |
818 |
mmeineke |
821 |
|
819 |
|
|
// create an fcc lattice in the water box. |
820 |
|
|
|
821 |
|
|
which = 0; |
822 |
|
|
for( i=0; i < nCellsX; i++ ){ |
823 |
|
|
for( j=0; j < nCellsY; j++ ){ |
824 |
|
|
for( k=0; k < nCellsZ; k++ ){ |
825 |
|
|
|
826 |
|
|
myORTHO->getLatticePoints(&posX, &posY, &posZ, i, j, k); |
827 |
|
|
for(l=0; l<4; l++){ |
828 |
|
|
waterSites[which].pos[0] = posX[l]; |
829 |
|
|
waterSites[which].pos[1] = posY[l]; |
830 |
|
|
waterSites[which].pos[2] = posZ[l]; |
831 |
mmeineke |
825 |
getRandomRot( waterSites[which].rot ); |
832 |
mmeineke |
821 |
which++; |
833 |
|
|
} |
834 |
|
|
} |
835 |
|
|
} |
836 |
|
|
} |
837 |
|
|
|
838 |
mmeineke |
825 |
// calc the system sizes |
839 |
mmeineke |
821 |
|
840 |
|
|
nAtoms = 0; |
841 |
|
|
molIndex = 0; |
842 |
mmeineke |
825 |
if(molStart != NULL ) delete[] molStart; |
843 |
mmeineke |
821 |
molStart = new int[nLipids + nWaters]; |
844 |
|
|
|
845 |
|
|
for(j=0; j<nLipids; j++){ |
846 |
|
|
molStart[molIndex] = nAtoms; |
847 |
|
|
molIndex++; |
848 |
|
|
nAtoms += lipidNatoms; |
849 |
|
|
} |
850 |
|
|
|
851 |
|
|
for(j=0; j<nWaters; j++){ |
852 |
|
|
molStart[molIndex] = nAtoms; |
853 |
|
|
molIndex++; |
854 |
|
|
nAtoms += waterNatoms; |
855 |
|
|
} |
856 |
|
|
|
857 |
mmeineke |
825 |
testInfo->n_atoms = nAtoms; |
858 |
|
|
theConfig = testInfo->getConfiguration(); |
859 |
|
|
theConfig->destroyArrays(); |
860 |
|
|
theConfig->createArrays( nAtoms ); |
861 |
|
|
testInfo->atoms = new Atom*[nAtoms]; |
862 |
|
|
atoms = testInfo->atoms; |
863 |
|
|
|
864 |
|
|
molIndex = 0; |
865 |
|
|
for(i=0; i<nLipids; i++ ){ |
866 |
|
|
lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms, |
867 |
|
|
molStart[molIndex], theConfig ); |
868 |
|
|
molIndex++; |
869 |
|
|
} |
870 |
|
|
|
871 |
|
|
for(i=0; i<nWaters; i++){ |
872 |
|
|
|
873 |
|
|
getRandomRot( waterSites[i].rot ); |
874 |
|
|
waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms, |
875 |
|
|
molStart[molIndex], theConfig ); |
876 |
|
|
molIndex++; |
877 |
|
|
} |
878 |
|
|
|
879 |
|
|
atoms[0]->getPos( myPos ); |
880 |
|
|
maxZ = myPos[2]; |
881 |
|
|
minZ = myPos[2]; |
882 |
|
|
for(i=0;i<nAtoms;i++){ |
883 |
|
|
atoms[i]->getPos( myPos ); |
884 |
|
|
minZ = (minZ > myPos[2]) ? myPos[2] : minZ; |
885 |
|
|
maxZ = (maxZ < myPos[2]) ? myPos[2] : maxZ; |
886 |
|
|
} |
887 |
|
|
boxZ = (maxZ - (minZ - waterFudge / 2.0)); |
888 |
|
|
|
889 |
|
|
// create the real Atom arrays |
890 |
|
|
|
891 |
mmeineke |
874 |
delete[] (mainInfo->atoms); |
892 |
|
|
theConfig = mainInfo->getConfiguration(); |
893 |
mmeineke |
821 |
theConfig->createArrays( nAtoms ); |
894 |
mmeineke |
874 |
mainInfo->atoms = new Atom*[nAtoms]; |
895 |
|
|
mainInfo->n_atoms = nAtoms; |
896 |
|
|
atoms = mainInfo->atoms; |
897 |
mmeineke |
821 |
|
898 |
|
|
// wrap back to <0,0,0> as center |
899 |
|
|
|
900 |
|
|
double Hmat[3][3]; |
901 |
|
|
|
902 |
|
|
Hmat[0][0] = boxX; |
903 |
|
|
Hmat[0][1] = 0.0; |
904 |
|
|
Hmat[0][2] = 0.0; |
905 |
|
|
|
906 |
|
|
Hmat[1][0] = 0.0; |
907 |
|
|
Hmat[1][1] = boxY; |
908 |
|
|
Hmat[1][2] = 0.0; |
909 |
|
|
|
910 |
|
|
Hmat[2][0] = 0.0; |
911 |
|
|
Hmat[2][1] = 0.0; |
912 |
|
|
Hmat[2][2] = boxZ; |
913 |
|
|
|
914 |
mmeineke |
825 |
mainInfo->setBoxM( Hmat ); |
915 |
mmeineke |
821 |
|
916 |
mmeineke |
825 |
for(j=0;j<nLipids;j++){ |
917 |
|
|
|
918 |
|
|
lipidSites[j].pos[0] -= centerX; |
919 |
|
|
lipidSites[j].pos[1] -= centerY; |
920 |
|
|
lipidSites[j].pos[2] -= centerZ; |
921 |
|
|
|
922 |
mmeineke |
874 |
mainInfo->wrapVector( lipidSites[j].pos ); |
923 |
mmeineke |
825 |
} |
924 |
mmeineke |
821 |
|
925 |
mmeineke |
825 |
for(j=0;j<nWaters;j++){ |
926 |
|
|
|
927 |
|
|
waterSites[j].pos[0] -= centerX; |
928 |
|
|
waterSites[j].pos[1] -= centerY; |
929 |
|
|
waterSites[j].pos[2] -= centerZ; |
930 |
|
|
|
931 |
mmeineke |
874 |
mainInfo->wrapVector( waterSites[j].pos ); |
932 |
mmeineke |
825 |
} |
933 |
mmeineke |
821 |
|
934 |
|
|
// initialize lipid positions |
935 |
|
|
|
936 |
|
|
molIndex = 0; |
937 |
|
|
for(i=0; i<nLipids; i++ ){ |
938 |
|
|
lipidLocate->placeMol( lipidSites[i].pos, lipidSites[i].rot, atoms, |
939 |
|
|
molStart[molIndex], theConfig ); |
940 |
|
|
molIndex++; |
941 |
|
|
} |
942 |
|
|
|
943 |
|
|
// initialize the water positions |
944 |
|
|
|
945 |
|
|
for(i=0; i<nWaters; i++){ |
946 |
mmeineke |
825 |
|
947 |
mmeineke |
821 |
getRandomRot( waterSites[i].rot ); |
948 |
|
|
waterLocate->placeMol( waterSites[i].pos, waterSites[i].rot, atoms, |
949 |
|
|
molStart[molIndex], theConfig ); |
950 |
|
|
molIndex++; |
951 |
|
|
} |
952 |
|
|
|
953 |
|
|
// set up the writer and write out |
954 |
|
|
|
955 |
mmeineke |
874 |
writer = new DumpWriter( mainInfo ); |
956 |
mmeineke |
821 |
writer->writeFinal( 0.0 ); |
957 |
|
|
|
958 |
mmeineke |
851 |
std::cout << "\n" |
959 |
|
|
<< "----------------------------------------\n" |
960 |
|
|
<< "\n" |
961 |
|
|
<< " final nLipids = " << nLipids << "\n" |
962 |
|
|
<< " final nWaters = " << nWaters << "\n" |
963 |
|
|
<< "\n"; |
964 |
|
|
|
965 |
mmeineke |
821 |
return 1; |
966 |
|
|
} |
967 |
|
|
|
968 |
|
|
|
969 |
|
|
/*************************************************************************** |
970 |
|
|
* prints out the usage for the command line arguments, then exits. |
971 |
|
|
***************************************************************************/ |
972 |
|
|
|
973 |
|
|
void usage(){ |
974 |
|
|
(void)fprintf(stdout, |
975 |
|
|
"\n" |
976 |
|
|
"The proper usage is: %s [options] <input_file>\n" |
977 |
|
|
"\n" |
978 |
|
|
"Options:\n" |
979 |
|
|
"\n" |
980 |
|
|
" short:\n" |
981 |
|
|
" ------\n" |
982 |
|
|
" -o <name> The output prefix\n" |
983 |
|
|
" -l <lipidName> The name of the lipid molecule specified in the BASS file.\n" |
984 |
|
|
" -w <waterName> The name of the water molecule specified in the BASS file.\n" |
985 |
|
|
" -s <double> The leaf spaceing of the bilayer.\n" |
986 |
|
|
" -a <double> The hexagonal lattice constant, a, for the bilayer leaf.\n" |
987 |
|
|
" -b <double> The hexagonal lattice constant, b, for the bilayer leaf.\n" |
988 |
|
|
" -h <double> Use to set a and b for a normal hex lattice with spacing <double>\n" |
989 |
|
|
|
990 |
|
|
"\n" |
991 |
|
|
" long:\n" |
992 |
|
|
" -----\n" |
993 |
|
|
|
994 |
|
|
" --version displays the version number\n" |
995 |
|
|
" --help displays this help message.\n" |
996 |
|
|
"\n" |
997 |
|
|
"\n", |
998 |
|
|
programName); |
999 |
|
|
} |