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